Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:10:55

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   double theElectronHcalIsoPt = 1.e+3;
0252   if (theElectron) {
0253     theElectronTrackIsoPt = theElectron->dr03TkSumPt();
0254     theElectronEcalIsoPt = theElectron->dr03EcalRecHitSumEt();
0255     theElectronHcalIsoPt = theElectron->dr03HcalTowerSumEt();
0256 
0257     if (electronIsoMode_ == kRelativeIso && theElectron->pt() > 0.) {
0258       theElectronTrackIsoPt /= theElectron->pt();
0259       theElectronEcalIsoPt /= theElectron->pt();
0260       theElectronHcalIsoPt /= theElectron->pt();
0261     }
0262   }
0263 
0264   //--- get collections of reconstructed tau-jets from the event

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

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

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

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

0384   // compute EWK tau analysis specific quantities

0385   //-----------------------------------------------------------------------------

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

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

0393   double mtElecPFMEt = calcMt(theElectron->px(), theElectron->py(), pfMEt.px(), pfMEt.py());
0394 
0395   // double pZetaCaloMEt = calcPzeta(theElectron->p4(), theTauJet->p4(),

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

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

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

0399 
0400   //-----------------------------------------------------------------------------

0401   // apply selection criteria; fill histograms

0402   //-----------------------------------------------------------------------------

0403 
0404   //--- fill electron multiplicity histogram

0405   unsigned numIdElectrons = 0;
0406   for (reco::GsfElectronCollection::const_iterator electron = electrons->begin(); electron != electrons->end();
0407        ++electron) {
0408     if (passesElectronId(*electron)) {
0409       ++numIdElectrons;
0410     }
0411   }
0412 
0413   // hNumIdElectrons_->Fill(numIdElectrons);

0414 
0415   ++numEventsAnalyzed_;
0416 
0417   bool isSelected = false;
0418   bool fullSelect = false;
0419   int cutFlowStatus = -1;
0420 
0421   if (mElecTau > visMassCut_) {
0422     cutFlowStatus = kPassedPreselection;
0423   }
0424   if (cutFlowStatus == kPassedPreselection && (isTriggered || hltPaths_.empty())) {
0425     cutFlowStatus = kPassedTrigger;
0426   }
0427   if (cutFlowStatus == kPassedTrigger && passesElectronId(*theElectron)) {
0428     cutFlowStatus = kPassedElectronId;
0429     hElectronTrackIsoPt_->Fill(theElectronTrackIsoPt);
0430   }
0431   if (cutFlowStatus == kPassedElectronId && theElectronTrackIsoPt < electronTrackIsoCut_) {
0432     cutFlowStatus = kPassedElectronTrackIso;
0433     hElectronEcalIsoPt_->Fill(theElectronEcalIsoPt);
0434   }
0435   if (cutFlowStatus == kPassedElectronTrackIso && theElectronEcalIsoPt < electronEcalIsoCut_) {
0436     cutFlowStatus = kPassedElectronEcalIso;
0437   }
0438   if (cutFlowStatus == kPassedElectronEcalIso && theTauDiscrByLeadTrackFinding > 0.5) {
0439     cutFlowStatus = kPassedTauLeadTrack;
0440     // if ( theTauJet->leadTrack().isAvailable() )

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

0442   }
0443   if (cutFlowStatus == kPassedTauLeadTrack && theTauDiscrByLeadTrackPtCut > 0.5) {
0444     cutFlowStatus = kPassedTauLeadTrackPt;
0445     // hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());

0446   }
0447   if (cutFlowStatus == kPassedTauLeadTrackPt && theTauDiscrAgainstElectrons > 0.5) {
0448     cutFlowStatus = kPassedTauDiscrAgainstElectrons;
0449     // hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);

0450   }
0451   if (cutFlowStatus == kPassedTauDiscrAgainstElectrons && theTauDiscrAgainstMuons > 0.5) {
0452     cutFlowStatus = kPassedTauDiscrAgainstMuons;
0453     isSelected = true;
0454   }
0455   if (cutFlowStatus == kPassedTauDiscrAgainstMuons && theTauDiscrByTrackIso > 0.5) {
0456     cutFlowStatus = kPassedTauTrackIso;
0457     // hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());

0458   }
0459   if (cutFlowStatus == kPassedTauTrackIso && theTauDiscrByEcalIso > 0.5) {
0460     cutFlowStatus = kPassedTauEcalIso;
0461     fullSelect = true;
0462     // hTauDiscrAgainstElectrons_->Fill(theTauDiscrAgainstElectrons);

0463   }
0464 
0465   for (int iCut = 1; iCut <= cutFlowStatus; ++iCut) {
0466     hCutFlowSummary_->Fill(iCut);
0467   }
0468 
0469   if (isSelected) {
0470     hElectronPt_->Fill(theElectron->pt());
0471     hElectronEta_->Fill(theElectron->eta());
0472     hElectronPhi_->Fill(theElectron->phi());
0473 
0474     hTauJetPt_->Fill(theTauJet->pt());
0475     hTauJetEta_->Fill(theTauJet->eta());
0476     // hTauJetPhi_->Fill(theTauJet->phi());

0477 
0478     // hTauJetCharge_->Fill(theTauJet->charge());

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

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

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

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

0483 
0484     if (fullSelect) {
0485       hVisMass_->Fill(mElecTau);
0486     }
0487     // hMtElecCaloMEt_->Fill(mtElecCaloMEt);

0488     hMtElecPFMEt_->Fill(mtElecPFMEt);
0489     // hPzetaCaloMEt_->Fill(pZetaCaloMEt);

0490     // hPzetaPFMEt_->Fill(pZetaPFMEt);

0491     hElecTauAcoplanarity_->Fill(dPhiElecTau);
0492     hElecTauCharge_->Fill(theElectron->charge() * theTauJet->charge());
0493 
0494     if (theEventVertex) {
0495       // hVertexChi2_->Fill(theEventVertex->normalizedChi2());

0496       hVertexZ_->Fill(theEventVertex->z());
0497       // hVertexD0_->Fill(getVertexD0(*theEventVertex, *beamSpot));

0498     }
0499 
0500     hCaloMEtPt_->Fill(caloMEt.pt());
0501     // hCaloMEtPhi_->Fill(caloMEt.phi());

0502 
0503     hPFMEtPt_->Fill(pfMEt.pt());
0504     // hPFMEtPhi_->Fill(pfMEt.phi());

0505   }
0506 
0507   if (isSelected)
0508     ++numEventsSelected_;
0509 }
0510 
0511 void EwkElecTauHistManager::finalizeHistograms() {
0512   edm::LogInfo("EwkElecTauHistManager") << "Filter-Statistics Summary:" << std::endl
0513                                         << " Events analyzed = " << numEventsAnalyzed_ << std::endl
0514                                         << " Events selected = " << numEventsSelected_;
0515   if (numEventsAnalyzed_ > 0) {
0516     double eff = numEventsSelected_ / (double)numEventsAnalyzed_;
0517     edm::LogInfo("") << "Overall efficiency = " << std::setprecision(4) << eff * 100. << " +/- " << std::setprecision(4)
0518                      << TMath::Sqrt(eff * (1 - eff) / numEventsAnalyzed_) * 100. << ")%";
0519   }
0520 }
0521 
0522 //-------------------------------------------------------------------------------

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

0524 //-------------------------------------------------------------------------------

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

0652   // access event-level information

0653   //-----------------------------------------------------------------------------

0654 
0655   bool readError = false;
0656 
0657   //--- get decision of high-level trigger for the event

0658   edm::Handle<edm::TriggerResults> hltDecision;
0659   readEventData(evt,
0660                 triggerResultsSource_,
0661                 hltDecision,
0662                 numWarningsTriggerResults_,
0663                 maxNumWarnings_,
0664                 readError,
0665                 "Failed to access Trigger results");
0666   if (readError)
0667     return;
0668 
0669   const edm::TriggerNames& triggerNames = evt.triggerNames(*hltDecision);
0670 
0671   bool isTriggered = false;
0672   for (vstring::const_iterator hltPath = hltPaths_.begin(); hltPath != hltPaths_.end(); ++hltPath) {
0673     unsigned int index = triggerNames.triggerIndex(*hltPath);
0674     if (index < triggerNames.size()) {
0675       if (hltDecision->accept(index))
0676         isTriggered = true;
0677     } else {
0678       if (numWarningsHLTpath_ < maxNumWarnings_ || maxNumWarnings_ == -1)
0679         edm::LogWarning("EwkMuTauHistManager") << " Undefined HLT path = " << (*hltPath) << " !!";
0680       ++numWarningsHLTpath_;
0681       continue;
0682     }
0683   }
0684 
0685   //--- get reconstructed primary event vertex of the event

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

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

0688   // Pt sum of tracks)

0689   edm::Handle<reco::VertexCollection> vertexCollection;
0690   readEventData(evt,
0691                 vertexSource_,
0692                 vertexCollection,
0693                 numWarningsVertex_,
0694                 maxNumWarnings_,
0695                 readError,
0696                 "Failed to access Vertex collection");
0697   if (readError)
0698     return;
0699 
0700   const reco::Vertex* theEventVertex = (!vertexCollection->empty()) ? &(vertexCollection->at(0)) : nullptr;
0701 
0702   //--- get beam-spot (expected vertex position) for the event

0703   edm::Handle<reco::BeamSpot> beamSpot;
0704   readEventData(
0705       evt, beamSpotSource_, beamSpot, numWarningsBeamSpot_, maxNumWarnings_, readError, "Failed to access Beam-spot");
0706   if (readError)
0707     return;
0708 
0709   //--- get collections of reconstructed muons from the event

0710   edm::Handle<reco::MuonCollection> muons;
0711   readEventData(
0712       evt, muonSource_, muons, numWarningsMuon_, maxNumWarnings_, readError, "Failed to access Muon collection");
0713   if (readError)
0714     return;
0715 
0716   const reco::Muon* theMuon = getTheMuon(*muons, muonEtaCut_, muonPtCut_);
0717 
0718   double theMuonTrackIsoPt = 1.e+3;
0719   double theMuonEcalIsoPt = 1.e+3;
0720   double theMuonCombIsoPt = 1.e+3;
0721 
0722   if (theMuon) {
0723     theMuonTrackIsoPt = theMuon->isolationR05().sumPt;
0724     // mu.isolationR05().emEt + mu.isolationR05().hadEt +

0725     // mu.isolationR05().sumPt

0726     theMuonEcalIsoPt = theMuon->isolationR05().emEt;
0727 
0728     if (muonIsoMode_ == kRelativeIso && theMuon->pt() > 0.) {
0729       theMuonTrackIsoPt /= theMuon->pt();
0730       theMuonEcalIsoPt /= theMuon->pt();
0731       theMuonCombIsoPt = (theMuon->isolationR05().sumPt + theMuon->isolationR05().emEt) / theMuon->pt();
0732       // std::cout<<"Rel Iso ="<<theMuonCombIsoPt<<std::endl;

0733     }
0734   }
0735 
0736   //--- get collections of reconstructed tau-jets from the event

0737   edm::Handle<reco::PFTauCollection> tauJets;
0738   readEventData(evt,
0739                 tauJetSource_,
0740                 tauJets,
0741                 numWarningsTauJet_,
0742                 maxNumWarnings_,
0743                 readError,
0744                 "Failed to access Tau-jet collection");
0745   if (readError)
0746     return;
0747 
0748   //--- get collections of tau-jet discriminators for those tau-jets

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

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

0815   edm::Handle<reco::CaloMETCollection> caloMEtCollection;
0816   readEventData(evt,
0817                 caloMEtSource_,
0818                 caloMEtCollection,
0819                 numWarningsCaloMEt_,
0820                 maxNumWarnings_,
0821                 readError,
0822                 "Failed to access calo. MET collection");
0823   if (readError)
0824     return;
0825 
0826   const reco::CaloMET& caloMEt = caloMEtCollection->at(0);
0827 
0828   edm::Handle<reco::PFMETCollection> pfMEtCollection;
0829   readEventData(evt,
0830                 pfMEtSource_,
0831                 pfMEtCollection,
0832                 numWarningsPFMEt_,
0833                 maxNumWarnings_,
0834                 readError,
0835                 "Failed to access pf. MET collection");
0836   if (readError)
0837     return;
0838 
0839   const reco::PFMET& pfMEt = pfMEtCollection->at(0);
0840 
0841   if (!(theMuon && theTauJet && theTauJetIndex != -1))
0842     return;
0843 
0844   //-----------------------------------------------------------------------------

0845   // compute EWK tau analysis specific quantities

0846   //-----------------------------------------------------------------------------

0847 
0848   double dPhiMuTau = calcDeltaPhi(theMuon->phi(), theTauJet->phi());
0849   // double dRMuTau = calcDeltaR(theMuon->p4(), theTauJet->p4());

0850   double dRMuTau = fabs(ROOT::Math::VectorUtil::DeltaR(theMuon->p4(), theTauJet->p4()));
0851   double mMuTau = (theMuon->p4() + theTauJet->p4()).M();
0852 
0853   // double mtMuCaloMEt = calcMt(theMuon->px(), theMuon->px(), caloMEt.px(),

0854   // caloMEt.py());

0855   double mtMuPFMEt = calcMt(theMuon->px(), theMuon->px(), pfMEt.px(), pfMEt.py());
0856 
0857   // double pZetaCaloMEt = calcPzeta(theMuon->p4(), theTauJet->p4(),

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

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

0860   // pfMEt.py());

0861 
0862   //-----------------------------------------------------------------------------

0863   // apply selection criteria; fill histograms

0864   //-----------------------------------------------------------------------------

0865 
0866   //--- fill muon multiplicity histogram

0867   unsigned numGlobalMuons = 0;
0868   for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0869     if (muon->isGlobalMuon()) {
0870       ++numGlobalMuons;
0871     }
0872   }
0873 
0874   // hNumGlobalMuons_->Fill(numGlobalMuons);

0875 
0876   ++numEventsAnalyzed_;
0877 
0878   bool isSelected = false;
0879   int cutFlowStatus = -1;
0880 
0881   // if ( muonIsoMode_ == kAbsoluteIso){

0882   if (mMuTau > visMassCut_) {
0883     cutFlowStatus = kPassedPreselection;
0884   }
0885   if (cutFlowStatus == kPassedPreselection && (isTriggered || hltPaths_.empty())) {
0886     cutFlowStatus = kPassedTrigger;
0887   }
0888   if (cutFlowStatus == kPassedTrigger && (theMuon->isGlobalMuon() || theMuon->isTrackerMuon())) {
0889     cutFlowStatus = kPassedMuonId;
0890   }
0891 
0892   if (cutFlowStatus == kPassedMuonId && (theTauDiscrByLeadTrackFinding > 0.5) && (theTauJet->eta() < tauJetEtaCut_) &&
0893       (theTauJet->pt() > tauJetPtCut_)) {
0894     cutFlowStatus = kPassedTauLeadTrack;
0895   }
0896   if (cutFlowStatus == kPassedTauLeadTrack && theTauDiscrByLeadTrackPtCut > 0.5) {
0897     cutFlowStatus = kPassedTauLeadTrackPt;
0898     // hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());

0899   }
0900   if (cutFlowStatus == kPassedTauLeadTrackPt && theTauDiscrAgainstMuons > 0.5) {
0901     cutFlowStatus = kPassedTauDiscrAgainstMuons;
0902     // hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());

0903   }
0904   if (cutFlowStatus == kPassedTauDiscrAgainstMuons && dRMuTau > deltaRCut_) {
0905     cutFlowStatus = kPassedDeltaR;
0906     // hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);

0907 
0908     hMuonPt_->Fill(theMuon->pt());
0909     hMuonEta_->Fill(theMuon->eta());
0910     hMuonPhi_->Fill(theMuon->phi());
0911 
0912     hTauJetPt_->Fill(theTauJet->pt());
0913     hTauJetEta_->Fill(theTauJet->eta());
0914     hTauJetPhi_->Fill(theTauJet->phi());
0915 
0916     // hTauJetCharge_->Fill(theTauJet->charge());

0917     if (theTauJet->signalTracks().isAvailable())
0918       hTauJetNumSignalTracks_->Fill(theTauJet->signalTracks().size());
0919     if (theTauJet->isolationTracks().isAvailable())
0920       hTauJetNumIsoTracks_->Fill(theTauJet->isolationTracks().size());
0921 
0922     hVisMass_->Fill(mMuTau);
0923     // hMtMuCaloMEt_->Fill(mtMuCaloMEt);

0924     hMtMuPFMEt_->Fill(mtMuPFMEt);
0925     // hPzetaCaloMEt_->Fill(pZetaCaloMEt);

0926     // hPzetaPFMEt_->Fill(pZetaPFMEt);

0927     hMuTauAcoplanarity_->Fill(dPhiMuTau);
0928     hMuTauDeltaR_->Fill(dRMuTau);
0929     // hMuTauCharge_->Fill(theMuon->charge() + theTauJet->charge());

0930 
0931     if (theEventVertex) {
0932       // hVertexChi2_->Fill(theEventVertex->normalizedChi2());

0933       hVertexZ_->Fill(theEventVertex->z());
0934       // hVertexD0_->Fill(getVertexD0(*theEventVertex, *beamSpot));

0935     }
0936 
0937     hCaloMEtPt_->Fill(caloMEt.pt());
0938     // hCaloMEtPhi_->Fill(caloMEt.phi());

0939 
0940     hPFMEtPt_->Fill(pfMEt.pt());
0941     // hPFMEtPhi_->Fill(pfMEt.phi());

0942     hMuonTrackIsoPt_->Fill(theMuonTrackIsoPt);
0943     hMuonEcalIsoPt_->Fill(theMuonEcalIsoPt);
0944     hMuonCombIsoPt_->Fill(theMuonCombIsoPt);
0945     // hMuonCombIsoPt_->Fill((theMuonTrackIsoPt+theMuonEcalIsoPt)/theMuon->pt());

0946 
0947     // std::cout<<"Rel Iso Hist =

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

0949     hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());
0950     hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());
0951     hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);
0952     if (theTauJet->leadTrack().isAvailable())
0953       hTauLeadTrackPt_->Fill(theTauJet->leadTrack()->pt());
0954   }
0955 
0956   if ((cutFlowStatus == kPassedDeltaR) && (((theMuonTrackIsoPt < muonTrackIsoCut_) && (muonIsoMode_ == kAbsoluteIso)) ||
0957                                            ((1 > 0) && (muonIsoMode_ == kRelativeIso)))) {
0958     cutFlowStatus = kPassedMuonTrackIso;
0959     // isSelected = true;

0960   }
0961   if (cutFlowStatus == kPassedMuonTrackIso &&
0962       (((theMuonEcalIsoPt < muonEcalIsoCut_) && (muonIsoMode_ == kAbsoluteIso)) ||
0963        ((theMuonCombIsoPt < muonCombIsoCut_) && (muonIsoMode_ == kRelativeIso)))) {
0964     cutFlowStatus = kPassedMuonEcalIso;
0965     // isSelected = true;

0966   }
0967 
0968   if (cutFlowStatus == kPassedMuonEcalIso && theTauDiscrByTrackIso > 0.5) {
0969     cutFlowStatus = kPassedTauTrackIso;
0970   }
0971 
0972   if (cutFlowStatus == kPassedTauTrackIso && theTauDiscrByEcalIso > 0.5) {
0973     cutFlowStatus = kPassedTauEcalIso;
0974     isSelected = true;
0975   }
0976 
0977   for (int iCut = 1; iCut <= cutFlowStatus; ++iCut) {
0978     hCutFlowSummary_->Fill(iCut);
0979   }
0980 
0981   for (int iCut = 1; iCut <= cutFlowStatus; ++iCut) {
0982     hCutFlowSummary_->Fill(iCut);
0983   }
0984 
0985   //     }

0986 
0987   if (isSelected) {
0988     hVisMassFinal_->Fill(mMuTau);
0989     ++numEventsSelected_;
0990   }
0991 }
0992 
0993 void EwkMuTauHistManager::finalizeHistograms() {
0994   edm::LogInfo("EwkMuTauHistManager") << "Filter-Statistics Summary:" << std::endl
0995                                       << " Events analyzed = " << numEventsAnalyzed_ << std::endl
0996                                       << " Events selected = " << numEventsSelected_;
0997   if (numEventsAnalyzed_ > 0) {
0998     double eff = numEventsSelected_ / (double)numEventsAnalyzed_;
0999     edm::LogInfo("") << "Overall efficiency = " << std::setprecision(4) << eff * 100. << " +/- " << std::setprecision(4)
1000                      << TMath::Sqrt(eff * (1 - eff) / numEventsAnalyzed_) * 100. << ")%";
1001   }
1002 }
1003 
1004 //-------------------------------------------------------------------------------

1005 // common auxiliary functions used by different channels

1006 //-------------------------------------------------------------------------------

1007 
1008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
1009 
1010 #include <TMath.h>
1011 
1012 int getIsoMode(const std::string& isoMode_string, int& error) {
1013   int isoMode_int;
1014   if (isoMode_string == "absoluteIso") {
1015     isoMode_int = kAbsoluteIso;
1016   } else if (isoMode_string == "relativeIso") {
1017     isoMode_int = kRelativeIso;
1018   } else {
1019     edm::LogError("getIsoMode") << " Failed to decode isoMode string = " << isoMode_string << " !!";
1020     isoMode_int = kUndefinedIso;
1021     error = 1;
1022   }
1023   return isoMode_int;
1024 }
1025 
1026 double calcDeltaPhi(double phi1, double phi2) {
1027   double deltaPhi = phi1 - phi2;
1028 
1029   if (deltaPhi < 0.)
1030     deltaPhi = -deltaPhi;
1031 
1032   if (deltaPhi > TMath::Pi())
1033     deltaPhi = 2 * TMath::Pi() - deltaPhi;
1034 
1035   return deltaPhi;
1036 }
1037 
1038 double calcMt(double px1, double py1, double px2, double py2) {
1039   double pt1 = TMath::Sqrt(px1 * px1 + py1 * py1);
1040   double pt2 = TMath::Sqrt(px2 * px2 + py2 * py2);
1041 
1042   double p1Dotp2 = px1 * px2 + py1 * py2;
1043   double cosAlpha = p1Dotp2 / (pt1 * pt2);
1044 
1045   return TMath::Sqrt(2 * pt1 * pt2 * (1 - cosAlpha));
1046 }
1047 
1048 double calcPzeta(const reco::Candidate::LorentzVector& p1,
1049                  const reco::Candidate::LorentzVector& p2,
1050                  double pxMEt,
1051                  double pyMEt) {
1052   double cosPhi1 = cos(p1.phi());
1053   double sinPhi1 = sin(p1.phi());
1054   double cosPhi2 = cos(p2.phi());
1055   double sinPhi2 = sin(p2.phi());
1056   double zetaX = cosPhi1 + cosPhi2;
1057   double zetaY = sinPhi1 + sinPhi2;
1058   double zetaR = TMath::Sqrt(zetaX * zetaX + zetaY * zetaY);
1059   if (zetaR > 0.) {
1060     zetaX /= zetaR;
1061     zetaY /= zetaR;
1062   }
1063 
1064   double pxVis = p1.px() + p2.px();
1065   double pyVis = p1.py() + p2.py();
1066   double pZetaVis = pxVis * zetaX + pyVis * zetaY;
1067 
1068   double px = pxVis + pxMEt;
1069   double py = pyVis + pyMEt;
1070   double pZeta = px * zetaX + py * zetaY;
1071 
1072   return pZeta - 1.5 * pZetaVis;
1073 }
1074 
1075 bool passesElectronPreId(const reco::GsfElectron& electron) {
1076   if ((TMath::Abs(electron.eta()) < 1.479 || TMath::Abs(electron.eta()) > 1.566) &&  // cut ECAL barrel/endcap crack

1077       electron.deltaPhiSuperClusterTrackAtVtx() < 0.8 && electron.deltaEtaSuperClusterTrackAtVtx() < 0.01 &&
1078       electron.sigmaIetaIeta() < 0.03) {
1079     return true;
1080   } else {
1081     return false;
1082   }
1083 }
1084 
1085 bool passesElectronId(const reco::GsfElectron& electron) {
1086   if (passesElectronPreId(electron) && ((TMath::Abs(electron.eta()) > 1.566 &&  // electron reconstructed in ECAL

1087                                                                                 // endcap

1088                                          electron.sigmaEtaEta() < 0.03 && electron.hcalOverEcal() < 0.05 &&
1089                                          TMath::Abs(electron.deltaEtaSuperClusterTrackAtVtx()) < 0.009 &&
1090                                          TMath::Abs(electron.deltaPhiSuperClusterTrackAtVtx()) < 0.7) ||
1091                                         (TMath::Abs(electron.eta()) < 1.479 &&  // electron reconstructed in ECAL

1092                                                                                 // barrel

1093                                          electron.sigmaEtaEta() < 0.01 && electron.hcalOverEcal() < 0.12 &&
1094                                          TMath::Abs(electron.deltaEtaSuperClusterTrackAtVtx()) < 0.007 &&
1095                                          TMath::Abs(electron.deltaPhiSuperClusterTrackAtVtx()) < 0.8))) {
1096     return true;
1097   } else {
1098     return false;
1099   }
1100 }
1101 
1102 const reco::GsfElectron* getTheElectron(const reco::GsfElectronCollection& electrons,
1103                                         double electronEtaCut,
1104                                         double electronPtCut) {
1105   const reco::GsfElectron* theElectron = nullptr;
1106 
1107   for (reco::GsfElectronCollection::const_iterator electron = electrons.begin(); electron != electrons.end();
1108        ++electron) {
1109     if (TMath::Abs(electron->eta()) < electronEtaCut && electron->pt() > electronPtCut &&
1110         passesElectronPreId(*electron)) {
1111       if (theElectron == nullptr || electron->pt() > theElectron->pt())
1112         theElectron = &(*electron);
1113     }
1114   }
1115 
1116   return theElectron;
1117 }
1118 
1119 const reco::Muon* getTheMuon(const reco::MuonCollection& muons, double muonEtaCut, double muonPtCut) {
1120   const reco::Muon* theMuon = nullptr;
1121 
1122   for (reco::MuonCollection::const_iterator muon = muons.begin(); muon != muons.end(); ++muon) {
1123     if (TMath::Abs(muon->eta()) < muonEtaCut && muon->pt() > muonPtCut) {
1124       if (theMuon == nullptr || muon->pt() > theMuon->pt())
1125         theMuon = &(*muon);
1126     }
1127   }
1128 
1129   return theMuon;
1130 }
1131 
1132 const reco::PFTau* getTheTauJet(const reco::PFTauCollection& tauJets,
1133                                 double tauJetEtaCut,
1134                                 double tauJetPtCut,
1135                                 int& theTauJetIndex) {
1136   const reco::PFTau* theTauJet = nullptr;
1137   theTauJetIndex = -1;
1138 
1139   int numTauJets = tauJets.size();
1140   for (int iTauJet = 0; iTauJet < numTauJets; ++iTauJet) {
1141     const reco::PFTau& tauJet = tauJets.at(iTauJet);
1142 
1143     if (fabs(tauJet.eta()) < tauJetEtaCut && tauJet.pt() > tauJetPtCut) {
1144       if (theTauJet == nullptr || tauJet.pt() > theTauJet->pt()) {
1145         theTauJet = &tauJet;
1146         theTauJetIndex = iTauJet;
1147       }
1148     }
1149   }
1150 
1151   return theTauJet;
1152 }
1153 
1154 double getVertexD0(const reco::Vertex& vertex, const reco::BeamSpot& beamSpot) {
1155   double dX = vertex.x() - beamSpot.x0();
1156   double dY = vertex.y() - beamSpot.y0();
1157   return TMath::Sqrt(dX * dX + dY * dY);
1158 }