Reference Manual

Go to the documentation of this file.
00001 #include "DQM/Physics/src/EwkTauDQM.h"
00002 
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 
00005 #include "FWCore/ServiceRegistry/interface/Service.h"
00006 #include "DQMServices/Core/interface/DQMStore.h"
00007 
00008 const std::string dqmSeparator = "/";
00009 
00010 std::string dqmDirectoryName(const std::string& dqmRootDirectory, const std::string& dqmSubDirectory)
00011 {
00012 //--- concatenate names of dqmRootDirectory and dqmSubDirectory;
00013 //    add "/" separator inbetween if necessary
00014   std::string dirName = dqmRootDirectory;
00015   if ( dirName != "" && dirName.find_last_of(dqmSeparator) != (dirName.length() - 1) )  dirName.append(dqmSeparator);
00016   dirName.append(dqmSubDirectory);
00017   return dirName;
00018 }
00019 
00020 EwkTauDQM::EwkTauDQM(const edm::ParameterSet& cfg)
00021   : dqmDirectory_(cfg.getParameter<std::string>("dqmDirectory")),
00022     dqmError_(0)
00023 {
00024   if ( !edm::Service<DQMStore>().isAvailable() ) {
00025     edm::LogError ("EwkTauDQM") << " Failed to access dqmStore --> histograms will NEITHER be booked NOR filled !!";
00026     dqmError_ = 1;
00027     return;
00028   }
00029 
00030   DQMStore* dqmStore = &(*edm::Service<DQMStore>());
00031 
00032   maxNumWarnings_ = cfg.exists("maxNumWarnings") ? cfg.getParameter<int>("maxNumWarnings") : 1;
00033 
00034   edm::ParameterSet cfgChannels = cfg.getParameter<edm::ParameterSet>("channels");
00035 
00036   edm::ParameterSet cfgElecTauChannel = cfgChannels.getParameter<edm::ParameterSet>("elecTauChannel");
00037   std::string dqmSubDirectoryElecTauChannel = cfgElecTauChannel.getParameter<std::string>("dqmSubDirectory");  
00038   cfgElecTauChannel.addParameter<std::string>("dqmDirectory", dqmDirectoryName(dqmDirectory_, dqmSubDirectoryElecTauChannel));
00039   cfgElecTauChannel.addParameter<int>("maxNumWarnings", maxNumWarnings_);
00040   elecTauHistManager_ = new EwkElecTauHistManager(cfgElecTauChannel, dqmStore);
00041 
00042   edm::ParameterSet cfgMuTauChannel = cfgChannels.getParameter<edm::ParameterSet>("muTauChannel");
00043   std::string dqmSubDirectoryMuTauChannel = cfgMuTauChannel.getParameter<std::string>("dqmSubDirectory");
00044   cfgMuTauChannel.addParameter<std::string>("dqmDirectory", dqmDirectoryName(dqmDirectory_, dqmSubDirectoryMuTauChannel));
00045   cfgMuTauChannel.addParameter<int>("maxNumWarnings", maxNumWarnings_);
00046   muTauHistManager_ = new EwkMuTauHistManager(cfgMuTauChannel, dqmStore);
00047 }
00048 
00049 EwkTauDQM::~EwkTauDQM()
00050 {
00051   delete elecTauHistManager_;
00052   delete muTauHistManager_;
00053 }
00054 
00055 void EwkTauDQM::beginJob()
00056 {
00057   if ( dqmError_ ) return;
00058  
00059   elecTauHistManager_->bookHistograms();
00060   muTauHistManager_->bookHistograms();
00061 }
00062 
00063 void EwkTauDQM::analyze(const edm::Event& evt, const edm::EventSetup& es)
00064 {
00065   if ( dqmError_ ) return;
00066 
00067   elecTauHistManager_->fillHistograms(evt, es);
00068   muTauHistManager_->fillHistograms(evt, es);
00069 }
00070 
00071 void EwkTauDQM::endJob()
00072 {
00073   if ( dqmError_ ) return;
00074 
00075   elecTauHistManager_->finalizeHistograms();
00076   muTauHistManager_->finalizeHistograms();
00077 }
00078 
00079 
00080 
00081 
00082 //-------------------------------------------------------------------------------
00083 // code specific to Z --> e + tau-jet channel
00084 //-------------------------------------------------------------------------------
00085 
00086 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00087 
00088 #include "FWCore/ServiceRegistry/interface/Service.h"
00089 #include "DQMServices/Core/interface/DQMStore.h"
00090 
00091 #include "DataFormats/Common/interface/Handle.h"
00092 #include "DataFormats/Common/interface/View.h"
00093 #include "DataFormats/Common/interface/TriggerResults.h"
00094 #include "FWCore/Common/interface/TriggerNames.h"
00095 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00096 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00097 #include "DataFormats/TauReco/interface/PFTau.h"
00098 #include "DataFormats/TauReco/interface/PFTauFwd.h"
00099 #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
00100 #include "DataFormats/METReco/interface/CaloMET.h"
00101 #include "DataFormats/METReco/interface/CaloMETFwd.h"
00102 #include "DataFormats/METReco/interface/PFMET.h"
00103 #include "DataFormats/METReco/interface/PFMETFwd.h"
00104 #include "DataFormats/TrackReco/interface/Track.h"
00105 #include "DataFormats/VertexReco/interface/Vertex.h"
00106 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00107 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00108 
00109 #include "TMath.h"
00110 
00111 #include <iostream>
00112 #include <iomanip>
00113 
00114 EwkElecTauHistManager::EwkElecTauHistManager(const edm::ParameterSet& cfg, DQMStore* dqmStore) 
00115   : dqmStore_(dqmStore),
00116     dqmDirectory_(cfg.getParameter<std::string>("dqmDirectory")),
00117     numEventsAnalyzed_(0),
00118     numEventsSelected_(0),
00119     cfgError_(0),
00120     numWarningsTriggerResults_(0),
00121     numWarningsHLTpath_(0),
00122     numWarningsVertex_(0),
00123     numWarningsBeamSpot_(0),
00124     numWarningsElectron_(0),
00125     numWarningsTauJet_(0),
00126     numWarningsTauDiscrByLeadTrackFinding_(0),
00127     numWarningsTauDiscrByLeadTrackPtCut_(0),
00128     numWarningsTauDiscrByTrackIso_(0),
00129     numWarningsTauDiscrByEcalIso_(0),
00130     numWarningsTauDiscrAgainstElectrons_(0),
00131     numWarningsTauDiscrAgainstMuons_(0),
00132     numWarningsCaloMEt_(0),
00133     numWarningsPFMEt_(0)
00134 {
00135   triggerResultsSource_ = cfg.getParameter<edm::InputTag>("triggerResultsSource");
00136   vertexSource_ = cfg.getParameter<edm::InputTag>("vertexSource");
00137   beamSpotSource_ = cfg.getParameter<edm::InputTag>("beamSpotSource");
00138   electronSource_ = cfg.getParameter<edm::InputTag>("electronSource");
00139   tauJetSource_ = cfg.getParameter<edm::InputTag>("tauJetSource");
00140   caloMEtSource_ = cfg.getParameter<edm::InputTag>("caloMEtSource");
00141   pfMEtSource_ = cfg.getParameter<edm::InputTag>("pfMEtSource");
00142 
00143   tauDiscrByLeadTrackFinding_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackFinding");
00144   tauDiscrByLeadTrackPtCut_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackPtCut");
00145   tauDiscrByTrackIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByTrackIso");
00146   tauDiscrByEcalIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByEcalIso");
00147   tauDiscrAgainstElectrons_ = cfg.getParameter<edm::InputTag>("tauDiscrAgainstElectrons");
00148   tauDiscrAgainstMuons_ = cfg.getParameter<edm::InputTag>("tauDiscrAgainstMuons");
00149 
00150   hltPaths_ = cfg.getParameter<vstring>("hltPaths");
00151 
00152   electronEtaCut_ = cfg.getParameter<double>("electronEtaCut");
00153   electronPtCut_ = cfg.getParameter<double>("electronPtCut");
00154   electronTrackIsoCut_ = cfg.getParameter<double>("electronTrackIsoCut");
00155   electronEcalIsoCut_ = cfg.getParameter<double>("electronEcalIsoCut");
00156   std::string electronIsoMode_string = cfg.getParameter<std::string>("electronIsoMode");
00157   electronIsoMode_ = getIsoMode(electronIsoMode_string, cfgError_);
00158 
00159   tauJetEtaCut_ = cfg.getParameter<double>("tauJetEtaCut");
00160   tauJetPtCut_ = cfg.getParameter<double>("tauJetPtCut");
00161 
00162   visMassCut_ = cfg.getParameter<double>("visMassCut");
00163 
00164   maxNumWarnings_ = cfg.exists("maxNumWarnings") ? cfg.getParameter<int>("maxNumWarnings") : 1;
00165 }
00166 
00167 void EwkElecTauHistManager::bookHistograms()
00168 {
00169   dqmStore_->setCurrentFolder(dqmDirectory_);
00170 
00171   //hNumIdElectrons_ = dqmStore_->book1D("NumIdElectronsMuons" , "Num. id. Muons", 5, -0.5, 4.5);
00172   hElectronPt_ = dqmStore_->book1D("ElectronPt" , "P_{T}^{e}", 20, 0., 100.);
00173   hElectronEta_ = dqmStore_->book1D("ElectronEta" , "#eta_{e}", 20, -4.0, +4.0);
00174   hElectronPhi_ = dqmStore_->book1D("ElectronPhi" , "#phi_{e}", 20, -TMath::Pi(), +TMath::Pi());
00175   hElectronTrackIsoPt_ = dqmStore_->book1D("ElectronTrackIsoPt" , "Electron Track Iso.", 20, -0.01, 0.5);
00176   hElectronEcalIsoPt_ = dqmStore_->book1D("ElectronEcalIsoPt" , "Electron Ecal Iso.", 20, -0.01, 0.5);
00177   //hElectronHcalIsoPt_ = dqmStore_->book1D("ElectronHcalIsoPt" , "Electron Hcal Iso.", 20, -0.01, 0.5);
00178 
00179   hTauJetPt_ = dqmStore_->book1D("TauJetPt" , "P_{T}^{#tau-Jet}", 20, 0., 100.);
00180   hTauJetEta_ = dqmStore_->book1D("TauJetEta" , "#eta_{#tau-Jet}", 20, -4.0, +4.0);
00181   //hTauJetPhi_ = dqmStore_->book1D("TauJetPhi" , "#phi_{#tau-Jet}", 20, -TMath::Pi(), +TMath::Pi());
00182   //hTauLeadTrackPt_ = dqmStore_->book1D("TauLeadTrackPt" , "P_{T}^{#tau-Jet}", 20, 0., 50.);
00183   //hTauTrackIsoPt_ = dqmStore_->book1D("TauTrackIsoPt" , "Tau Track Iso.", 20, -0.01, 40.);
00184   //hTauEcalIsoPt_ = dqmStore_->book1D("TauEcalIsoPt" , "Tau Ecal Iso.", 10, -0.01, 10.);
00185   //hTauDiscrAgainstElectrons_ = dqmStore_->book1D("TauDiscrAgainstElectrons" , "Tau Discr. against Electrons", 2, -0.5, +1.5);
00186   //hTauDiscrAgainstMuons_ = dqmStore_->book1D("TauDiscrAgainstMuons" , "Tau Discr. against Muons", 2, -0.5, +1.5);
00187   //hTauJetCharge_ = dqmStore_->book1D("TauJetCharge" , "Q_{#tau-Jet}", 11, -5.5, +5.5);
00188   //hTauJetNumSignalTracks_ = dqmStore_->book1D("TauJetNumSignalTracks" , "Num. Tau signal Cone Tracks", 20, -0.5, +19.5);
00189   //hTauJetNumIsoTracks_ = dqmStore_->book1D("TauJetNumIsoTracks" , "Num. Tau isolation Cone Tracks", 20, -0.5, +19.5);
00190   
00191   hVisMass_ = dqmStore_->book1D("VisMass", "e + #tau-Jet visible Mass", 20, 20., 120.);
00192   //hMtElecCaloMEt_ = dqmStore_->book1D("MtElecCaloMEt", "e + E_{T}^{miss} (Calo) transverse Mass", 20, 20., 120.);
00193   hMtElecPFMEt_ = dqmStore_->book1D("MtElecPFMEt", "e + E_{T}^{miss} (PF) transverse Mass", 20, 20., 120.);
00194   //hPzetaCaloMEt_ = dqmStore_->book1D("PzetaCaloMEt", "P_{#zeta} - 1.5*P_{#zeta}^{vis} (Calo)", 20, -40., 40.);
00195   //hPzetaPFMEt_ = dqmStore_->book1D("PzetaPFMEt", "P_{#zeta} - 1.5*P_{#zeta}^{vis} (PF)", 20, -40., 40.);
00196   hElecTauAcoplanarity_ = dqmStore_->book1D("ElecTauAcoplanarity", "#Delta #phi_{e #tau-Jet}", 20, -TMath::Pi(), +TMath::Pi());
00197   hElecTauCharge_ = dqmStore_->book1D("ElecTauCharge" , "Q_{e * #tau-Jet}", 5, -2.5, +2.5);
00198 
00199   //hVertexChi2_ = dqmStore_->book1D("VertexChi2", "Event Vertex #chi^{2} / n.d.o.f.", 20, 0., 2.0);
00200   hVertexZ_ = dqmStore_->book1D("VertexZ", "Event Vertex z-Position", 20, -25., +25.);
00201   //hVertexD0_ = dqmStore_->book1D("VertexD0", "Event Vertex d_{0}", 20, -0.0001, 0.05);
00202 
00203   hCaloMEtPt_ = dqmStore_->book1D("CaloMEtPt", "E_{T}^{miss} (Calo)", 20, 0., 100.);
00204   //hCaloMEtPhi_ = dqmStore_->book1D("CaloMEtPhi", "#phi^{miss} (Calo)", 20, -TMath::Pi(), +TMath::Pi());
00205 
00206   hPFMEtPt_ = dqmStore_->book1D("PFMEtPt", "E_{T}^{miss} (PF)", 20, 0., 100.);
00207   //hPFMEtPhi_ = dqmStore_->book1D("PFMEtPhi", "#phi^{miss} (PF)", 20, -TMath::Pi(), +TMath::Pi());
00208 
00209   hCutFlowSummary_ = dqmStore_->book1D("CutFlowSummary", "Cut-flow Summary", 11, 0.5, 11.5);
00210   hCutFlowSummary_->setBinLabel(kPassedPreselection, "Preselection");
00211   hCutFlowSummary_->setBinLabel(kPassedTrigger, "HLT");
00212   hCutFlowSummary_->setBinLabel(kPassedElectronId, "e ID");
00213   hCutFlowSummary_->setBinLabel(kPassedElectronTrackIso, "e Trk Iso.");
00214   hCutFlowSummary_->setBinLabel(kPassedElectronEcalIso, "e Ecal Iso.");
00215   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrack, "#tau lead. Track");
00216   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrackPt, "#tau lead. Track P_{T}");
00217   hCutFlowSummary_->setBinLabel(kPassedTauDiscrAgainstElectrons, "#tau anti-e Discr.");
00218   hCutFlowSummary_->setBinLabel(kPassedTauDiscrAgainstMuons, "#tau anti-#mu Discr.");
00219   hCutFlowSummary_->setBinLabel(kPassedTauTrackIso, "#tau Track Iso.");
00220   hCutFlowSummary_->setBinLabel(kPassedTauEcalIso, "#tau Ecal Iso.");
00221 }
00222 
00223 void EwkElecTauHistManager::fillHistograms(const edm::Event& evt, const edm::EventSetup& es)
00224 {
00225   if ( cfgError_ ) return;
00226 
00227   //-----------------------------------------------------------------------------
00228   // access event-level information
00229   //-----------------------------------------------------------------------------
00230 
00231   bool readError = false;
00232 
00233 //--- get decision of high-level trigger for the event
00234   edm::Handle<edm::TriggerResults> hltDecision;
00235   readEventData(evt, triggerResultsSource_, hltDecision, numWarningsTriggerResults_, maxNumWarnings_, 
00236                 readError, "Failed to access Trigger results");
00237   if ( readError ) return;
00238   
00239   const edm::TriggerNames & triggerNames = evt.triggerNames(*hltDecision);
00240    
00241   bool isTriggered = false;
00242   for ( vstring::const_iterator hltPath = hltPaths_.begin();
00243         hltPath != hltPaths_.end(); ++hltPath ) {
00244     unsigned int index = triggerNames.triggerIndex(*hltPath);
00245     if ( index < triggerNames.size() ) {
00246       if ( hltDecision->accept(index) ) isTriggered = true;
00247     } else {
00248       if ( numWarningsHLTpath_ < maxNumWarnings_ || maxNumWarnings_ == -1 ) 
00249         edm::LogWarning ("EwkElecTauHistManager") << " Undefined HLT path = " << (*hltPath) << " !!";
00250       ++numWarningsHLTpath_;
00251       continue;
00252     }
00253   }
00254   
00255 //--- get reconstructed primary event vertex of the event
00256 //   (take as "the" primary event vertex the first entry in the collection
00257 //    of vertex objects, corresponding to the vertex associated to the highest Pt sum of tracks)
00258   edm::Handle<reco::VertexCollection> vertexCollection;
00259   readEventData(evt, vertexSource_, vertexCollection, numWarningsVertex_, maxNumWarnings_,
00260                 readError, "Failed to access Vertex collection");
00261   if ( readError ) return;
00262 
00263   const reco::Vertex* theEventVertex = ( vertexCollection->size() > 0 ) ? &(vertexCollection->at(0)) : 0;
00264 
00265 //--- get beam-spot (expected vertex position) for the event
00266   edm::Handle<reco::BeamSpot> beamSpot;
00267   readEventData(evt, beamSpotSource_, beamSpot, numWarningsBeamSpot_, maxNumWarnings_,
00268                 readError, "Failed to access Beam-spot");
00269   if ( readError ) return;
00270   
00271 //--- get collections of reconstructed electrons from the event
00272   edm::Handle<reco::GsfElectronCollection> electrons;
00273   readEventData(evt, electronSource_, electrons, numWarningsElectron_, maxNumWarnings_,
00274                 readError, "Failed to access Electron collection");
00275   if ( readError ) return;
00276 
00277   const reco::GsfElectron* theElectron = getTheElectron(*electrons, electronEtaCut_, electronPtCut_);
00278 
00279   double theElectronTrackIsoPt = 1.e+3;
00280   double theElectronEcalIsoPt = 1.e+3;
00281   double theElectronHcalIsoPt = 1.e+3;
00282   if ( theElectron ) {
00283     theElectronTrackIsoPt = theElectron->dr03TkSumPt();
00284     theElectronEcalIsoPt = theElectron->dr03EcalRecHitSumEt();
00285     theElectronHcalIsoPt = theElectron->dr03HcalTowerSumEt();
00286 
00287     if ( electronIsoMode_ == kRelativeIso && theElectron->pt() > 0. ) {
00288       theElectronTrackIsoPt /= theElectron->pt();
00289       theElectronEcalIsoPt /= theElectron->pt();
00290       theElectronHcalIsoPt /= theElectron->pt();
00291     }
00292   }
00293 
00294 //--- get collections of reconstructed tau-jets from the event
00295   edm::Handle<reco::PFTauCollection> tauJets;
00296   readEventData(evt, tauJetSource_, tauJets, numWarningsTauJet_, maxNumWarnings_,
00297                 readError, "Failed to access Tau-jet collection");
00298   if ( readError ) return;
00299 
00300 //--- get collections of tau-jet discriminators for those tau-jets
00301   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackFinding;
00302   readEventData(evt, tauDiscrByLeadTrackFinding_, tauDiscrByLeadTrackFinding, numWarningsTauDiscrByLeadTrackFinding_, maxNumWarnings_,
00303                 readError, "Failed to access collection of pf. Tau discriminators by leading Track finding");
00304   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackPtCut;
00305   readEventData(evt, tauDiscrByLeadTrackPtCut_, tauDiscrByLeadTrackPtCut, numWarningsTauDiscrByLeadTrackPtCut_, maxNumWarnings_,
00306                 readError, "Failed to access collection of pf. Tau discriminators by leading Track Pt cut");
00307   edm::Handle<reco::PFTauDiscriminator> tauDiscrByTrackIso;
00308   readEventData(evt, tauDiscrByTrackIso_, tauDiscrByTrackIso, numWarningsTauDiscrByTrackIso_, maxNumWarnings_,
00309                 readError, "Failed to access collection of pf. Tau discriminators by Track isolation");
00310   edm::Handle<reco::PFTauDiscriminator> tauDiscrByEcalIso;
00311   readEventData(evt, tauDiscrByTrackIso_, tauDiscrByEcalIso, numWarningsTauDiscrByEcalIso_, maxNumWarnings_,
00312                 readError, "Failed to access collection of pf. Tau discriminators by ECAL isolation");
00313   edm::Handle<reco::PFTauDiscriminator> tauDiscrAgainstElectrons;
00314   readEventData(evt, tauDiscrAgainstElectrons_, tauDiscrAgainstElectrons, numWarningsTauDiscrAgainstElectrons_, maxNumWarnings_,
00315                 readError, "Failed to access collection of pf. Tau discriminators against Electrons");
00316   edm::Handle<reco::PFTauDiscriminator> tauDiscrAgainstMuons;
00317   readEventData(evt, tauDiscrAgainstMuons_, tauDiscrAgainstMuons, numWarningsTauDiscrAgainstMuons_, maxNumWarnings_,
00318                 readError, "Failed to access collection of pf. Tau discriminators against Muons");
00319   if ( readError ) return;
00320 
00321   int theTauJetIndex = -1;
00322   const reco::PFTau* theTauJet = getTheTauJet(*tauJets, tauJetEtaCut_, tauJetPtCut_, theTauJetIndex);
00323 
00324   double theTauDiscrByLeadTrackFinding = -1.;
00325   double theTauDiscrByLeadTrackPtCut = -1.;
00326   double theTauDiscrByTrackIso = -1.;
00327   double theTauDiscrByEcalIso = -1.;
00328   double theTauDiscrAgainstElectrons = -1.;
00329   double theTauDiscrAgainstMuons = -1.;
00330   if ( theTauJetIndex != -1 ) {
00331     reco::PFTauRef theTauJetRef(tauJets, theTauJetIndex);
00332     theTauDiscrByLeadTrackFinding = (*tauDiscrByLeadTrackFinding)[theTauJetRef];
00333     theTauDiscrByLeadTrackPtCut = (*tauDiscrByLeadTrackPtCut)[theTauJetRef];
00334     theTauDiscrByTrackIso = (*tauDiscrByTrackIso)[theTauJetRef];
00335     theTauDiscrByEcalIso = (*tauDiscrByEcalIso)[theTauJetRef];
00336     theTauDiscrAgainstElectrons = (*tauDiscrAgainstElectrons)[theTauJetRef];
00337     theTauDiscrAgainstMuons = (*tauDiscrAgainstMuons)[theTauJetRef];
00338   }
00339 
00340 //--- get missing transverse momentum
00341 //    measured by calorimeters/reconstructed by particle-flow algorithm
00342   edm::Handle<reco::CaloMETCollection> caloMEtCollection;
00343   readEventData(evt, caloMEtSource_, caloMEtCollection, numWarningsCaloMEt_, maxNumWarnings_,
00344                 readError, "Failed to access calo. MET collection");
00345   if ( readError ) return;
00346 
00347   const reco::CaloMET& caloMEt = caloMEtCollection->at(0);
00348   
00349   edm::Handle<reco::PFMETCollection> pfMEtCollection;
00350   readEventData(evt, pfMEtSource_, pfMEtCollection, numWarningsPFMEt_, maxNumWarnings_,
00351                 readError, "Failed to access pf. MET collection");
00352   if ( readError ) return;
00353 
00354   const reco::PFMET& pfMEt = pfMEtCollection->at(0);
00355 
00356   if ( !(theElectron && theTauJet && theTauJetIndex != -1) ) return;
00357 
00358   //-----------------------------------------------------------------------------
00359   // compute EWK tau analysis specific quantities
00360   //-----------------------------------------------------------------------------
00361 
00362   double dPhiElecTau = calcDeltaPhi(theElectron->phi(), theTauJet->phi());
00363 
00364   double mElecTau = (theElectron->p4() + theTauJet->p4()).M();
00365 
00366   //double mtElecCaloMEt = calcMt(theElectron->px(), theElectron->py(), caloMEt.px(), caloMEt.py());
00367   double mtElecPFMEt = calcMt(theElectron->px(), theElectron->py(), pfMEt.px(), pfMEt.py());
00368 
00369   //double pZetaCaloMEt = calcPzeta(theElectron->p4(), theTauJet->p4(), caloMEt.px(), caloMEt.py());
00370   //double pZetaPFMEt = calcPzeta(theElectron->p4(), theTauJet->p4(), pfMEt.px(), pfMEt.py());
00371 
00372   //-----------------------------------------------------------------------------
00373   // apply selection criteria; fill histograms
00374   //-----------------------------------------------------------------------------
00375 
00376 //--- fill electron multiplicity histogram
00377   unsigned numIdElectrons = 0;
00378   for ( reco::GsfElectronCollection::const_iterator electron = electrons->begin();
00379         electron != electrons->end(); ++electron ) {
00380     if ( passesElectronId(*electron) ) {
00381       ++numIdElectrons;
00382     }
00383   }
00384 
00385   //hNumIdElectrons_->Fill(numIdElectrons);
00386 
00387   ++numEventsAnalyzed_;
00388 
00389   bool isSelected = false;
00390   bool fullSelect = false;
00391   int cutFlowStatus = -1;
00392 
00393   if ( mElecTau > visMassCut_ ) {
00394     cutFlowStatus = kPassedPreselection;
00395   }
00396   if ( cutFlowStatus == kPassedPreselection && (isTriggered || hltPaths_.size() == 0) ) {
00397     cutFlowStatus = kPassedTrigger;
00398   }
00399   if ( cutFlowStatus == kPassedTrigger && passesElectronId(*theElectron) ) {
00400     cutFlowStatus = kPassedElectronId;
00401     hElectronTrackIsoPt_->Fill(theElectronTrackIsoPt);
00402   }
00403   if ( cutFlowStatus == kPassedElectronId && theElectronTrackIsoPt < electronTrackIsoCut_ ) {
00404     cutFlowStatus = kPassedElectronTrackIso;
00405     hElectronEcalIsoPt_->Fill(theElectronEcalIsoPt);
00406   }
00407   if ( cutFlowStatus == kPassedElectronTrackIso && theElectronEcalIsoPt < electronEcalIsoCut_ ) {
00408     cutFlowStatus = kPassedElectronEcalIso;
00409   }
00410   if ( cutFlowStatus == kPassedElectronEcalIso && theTauDiscrByLeadTrackFinding > 0.5 ) {
00411     cutFlowStatus = kPassedTauLeadTrack;
00412     //if ( theTauJet->leadTrack().isAvailable() ) hTauLeadTrackPt_->Fill(theTauJet->leadTrack()->pt());
00413   }
00414   if ( cutFlowStatus == kPassedTauLeadTrack && theTauDiscrByLeadTrackPtCut > 0.5 ) {
00415     cutFlowStatus = kPassedTauLeadTrackPt;
00416     //hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());
00417   }
00418   if ( cutFlowStatus == kPassedTauLeadTrackPt && theTauDiscrAgainstElectrons > 0.5 ) {
00419     cutFlowStatus = kPassedTauDiscrAgainstElectrons;
00420     //hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);
00421   }
00422   if ( cutFlowStatus == kPassedTauDiscrAgainstElectrons && theTauDiscrAgainstMuons > 0.5 ) {
00423     cutFlowStatus = kPassedTauDiscrAgainstMuons;
00424     isSelected = true;
00425   }
00426   if ( cutFlowStatus == kPassedTauDiscrAgainstMuons && theTauDiscrByTrackIso > 0.5 ) {
00427     cutFlowStatus = kPassedTauTrackIso;
00428     //hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());
00429   }
00430   if ( cutFlowStatus == kPassedTauTrackIso && theTauDiscrByEcalIso > 0.5 ) {
00431     cutFlowStatus = kPassedTauEcalIso;
00432     fullSelect = true;
00433     //hTauDiscrAgainstElectrons_->Fill(theTauDiscrAgainstElectrons);
00434   }
00435 
00436 
00437   for ( int iCut = 1; iCut <= cutFlowStatus; ++iCut ) {
00438     hCutFlowSummary_->Fill(iCut);
00439   }
00440 
00441   if ( isSelected ) {
00442     hElectronPt_->Fill(theElectron->pt());
00443     hElectronEta_->Fill(theElectron->eta());
00444     hElectronPhi_->Fill(theElectron->phi());
00445 
00446     hTauJetPt_->Fill(theTauJet->pt());
00447     hTauJetEta_->Fill(theTauJet->eta());
00448     //hTauJetPhi_->Fill(theTauJet->phi());
00449 
00450     //hTauJetCharge_->Fill(theTauJet->charge());
00451     //if ( theTauJet->signalTracks().isAvailable()    ) hTauJetNumSignalTracks_->Fill(theTauJet->signalTracks().size());
00452     //if ( theTauJet->isolationTracks().isAvailable() ) hTauJetNumIsoTracks_->Fill(theTauJet->isolationTracks().size());
00453   
00454     if (fullSelect ) {hVisMass_->Fill(mElecTau);}
00455     //hMtElecCaloMEt_->Fill(mtElecCaloMEt);
00456     hMtElecPFMEt_->Fill(mtElecPFMEt);
00457     //hPzetaCaloMEt_->Fill(pZetaCaloMEt);
00458     //hPzetaPFMEt_->Fill(pZetaPFMEt);
00459     hElecTauAcoplanarity_->Fill(dPhiElecTau);
00460     hElecTauCharge_->Fill(theElectron->charge() * theTauJet->charge());
00461 
00462     if ( theEventVertex ) {
00463       //hVertexChi2_->Fill(theEventVertex->normalizedChi2());
00464       hVertexZ_->Fill(theEventVertex->z());
00465       //hVertexD0_->Fill(getVertexD0(*theEventVertex, *beamSpot));
00466     }
00467     
00468     hCaloMEtPt_->Fill(caloMEt.pt());
00469     //hCaloMEtPhi_->Fill(caloMEt.phi());
00470 
00471     hPFMEtPt_->Fill(pfMEt.pt());
00472     //hPFMEtPhi_->Fill(pfMEt.phi());
00473   }
00474 
00475   if ( isSelected ) ++numEventsSelected_;
00476 }
00477 
00478 void EwkElecTauHistManager::finalizeHistograms()
00479 {
00480   edm::LogInfo ("EwkElecTauHistManager") 
00481     << "Filter-Statistics Summary:" << std::endl
00482     << " Events analyzed = " << numEventsAnalyzed_ << std::endl
00483     << " Events selected = " << numEventsSelected_;
00484   if ( numEventsAnalyzed_ > 0 ) {
00485     double eff = numEventsSelected_/(double)numEventsAnalyzed_;
00486     edm::LogInfo ("") 
00487       << "Overall efficiency = " << std::setprecision(4) << eff*100. 
00488       << " +/- " << std::setprecision(4) << TMath::Sqrt(eff*(1 - eff)/numEventsAnalyzed_)*100. << ")%";
00489   }
00490 }
00491 
00492 
00493 
00494 
00495 //-------------------------------------------------------------------------------
00496 // code specific to Z --> mu + tau-jet channel
00497 //-------------------------------------------------------------------------------
00498 
00499 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00500 
00501 #include "FWCore/ServiceRegistry/interface/Service.h"
00502 #include "DQMServices/Core/interface/DQMStore.h"
00503 
00504 #include "DataFormats/Common/interface/Handle.h"
00505 #include "DataFormats/Common/interface/View.h"
00506 #include "DataFormats/Common/interface/TriggerResults.h"
00507 #include "FWCore/Common/interface/TriggerNames.h"
00508 #include "DataFormats/MuonReco/interface/Muon.h"
00509 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00510 #include "DataFormats/TauReco/interface/PFTau.h"
00511 #include "DataFormats/TauReco/interface/PFTauFwd.h"
00512 #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
00513 #include "DataFormats/METReco/interface/CaloMET.h"
00514 #include "DataFormats/METReco/interface/CaloMETFwd.h"
00515 #include "DataFormats/METReco/interface/PFMET.h"
00516 #include "DataFormats/METReco/interface/PFMETFwd.h"
00517 #include "DataFormats/TrackReco/interface/Track.h"
00518 #include "DataFormats/VertexReco/interface/Vertex.h"
00519 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00520 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00521 
00522 #include "TMath.h"
00523 
00524 #include <iostream>
00525 #include <iomanip>
00526 
00527 EwkMuTauHistManager::EwkMuTauHistManager(const edm::ParameterSet& cfg, DQMStore* dqmStore) 
00528   : dqmStore_(dqmStore),
00529     dqmDirectory_(cfg.getParameter<std::string>("dqmDirectory")),
00530     numEventsAnalyzed_(0),
00531     numEventsSelected_(0),
00532     cfgError_(0),
00533     numWarningsTriggerResults_(0),
00534     numWarningsHLTpath_(0),
00535     numWarningsVertex_(0),
00536     numWarningsBeamSpot_(0),
00537     numWarningsMuon_(0),
00538     numWarningsTauJet_(0),
00539     numWarningsTauDiscrByLeadTrackFinding_(0),
00540     numWarningsTauDiscrByLeadTrackPtCut_(0),
00541     numWarningsTauDiscrByTrackIso_(0),
00542     numWarningsTauDiscrByEcalIso_(0),
00543     numWarningsTauDiscrAgainstMuons_(0),
00544     numWarningsCaloMEt_(0),
00545     numWarningsPFMEt_(0)
00546 {
00547   triggerResultsSource_ = cfg.getParameter<edm::InputTag>("triggerResultsSource");
00548   vertexSource_ = cfg.getParameter<edm::InputTag>("vertexSource");
00549   beamSpotSource_ = cfg.getParameter<edm::InputTag>("beamSpotSource");
00550   muonSource_ = cfg.getParameter<edm::InputTag>("muonSource");
00551   tauJetSource_ = cfg.getParameter<edm::InputTag>("tauJetSource");
00552   caloMEtSource_ = cfg.getParameter<edm::InputTag>("caloMEtSource");
00553   pfMEtSource_ = cfg.getParameter<edm::InputTag>("pfMEtSource");
00554 
00555   tauDiscrByLeadTrackFinding_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackFinding");
00556   tauDiscrByLeadTrackPtCut_ = cfg.getParameter<edm::InputTag>("tauDiscrByLeadTrackPtCut");
00557   tauDiscrByTrackIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByTrackIso");
00558   tauDiscrByEcalIso_ = cfg.getParameter<edm::InputTag>("tauDiscrByEcalIso");
00559   tauDiscrAgainstMuons_ = cfg.getParameter<edm::InputTag>("tauDiscrAgainstMuons");
00560 
00561   hltPaths_ = cfg.getParameter<vstring>("hltPaths");
00562 
00563   muonEtaCut_ = cfg.getParameter<double>("muonEtaCut");
00564   muonPtCut_ = cfg.getParameter<double>("muonPtCut");
00565   muonTrackIsoCut_ = cfg.getParameter<double>("muonTrackIsoCut");
00566   muonEcalIsoCut_ = cfg.getParameter<double>("muonEcalIsoCut");
00567   muonCombIsoCut_ = cfg.getParameter<double>("muonCombIsoCut");
00568   std::string muonIsoMode_string = cfg.getParameter<std::string>("muonIsoMode");
00569   muonIsoMode_ = getIsoMode(muonIsoMode_string, cfgError_);
00570 
00571   tauJetEtaCut_ = cfg.getParameter<double>("tauJetEtaCut");
00572   tauJetPtCut_ = cfg.getParameter<double>("tauJetPtCut");
00573 
00574   visMassCut_ = cfg.getParameter<double>("visMassCut");
00575 deltaRCut_ = cfg.getParameter<double>("deltaRCut");
00576 
00577   maxNumWarnings_ = cfg.exists("maxNumWarnings") ? cfg.getParameter<int>("maxNumWarnings") : 1;
00578 }
00579 
00580 void EwkMuTauHistManager::bookHistograms()
00581 {
00582   dqmStore_->setCurrentFolder(dqmDirectory_);
00583 
00584   //hNumGlobalMuons_ = dqmStore_->book1D("NumGlobalMuons" , "Num. global Muons", 5, -0.5, 4.5);
00585   hMuonPt_ = dqmStore_->book1D("MuonPt" , "P_{T}^{#mu}", 20, 0., 100.);
00586   hMuonEta_ = dqmStore_->book1D("MuonEta" , "#eta_{#mu}", 20, -4.0, +4.0);
00587   hMuonPhi_ = dqmStore_->book1D("MuonPhi" , "#phi_{#mu}", 20, -TMath::Pi(), +TMath::Pi());
00588   hMuonTrackIsoPt_ = dqmStore_->book1D("MuonTrackIsoPt" , "Muon Track Iso.", 20, -0.01, 10.);
00589   hMuonEcalIsoPt_ = dqmStore_->book1D("MuonEcalIsoPt" , "Muon Ecal Iso.", 20, -0.01, 10.);
00590   hMuonCombIsoPt_ = dqmStore_->book1D("MuonCombIsoPt" , "Muon Comb Iso.", 20, -0.01, 1.);
00591 
00592   hTauJetPt_ = dqmStore_->book1D("TauJetPt" , "P_{T}^{#tau-Jet}", 20, 0., 100.);
00593   hTauJetEta_ = dqmStore_->book1D("TauJetEta" , "#eta_{#tau-Jet}", 20, -4.0, +4.0);
00594   hTauJetPhi_ = dqmStore_->book1D("TauJetPhi" , "#phi_{#tau-Jet}", 20, -TMath::Pi(), +TMath::Pi());
00595   hTauLeadTrackPt_ = dqmStore_->book1D("TauLeadTrackPt" , "P_{T}^{#tau-Jetldg trk}", 20, 0., 50.);
00596   hTauTrackIsoPt_ = dqmStore_->book1D("TauTrackIsoPt" , "Tau Track Iso.", 20, -0.01, 40.);
00597   hTauEcalIsoPt_ = dqmStore_->book1D("TauEcalIsoPt" , "Tau Ecal Iso.", 10, -0.01, 10.);
00598   hTauDiscrAgainstMuons_ = dqmStore_->book1D("TauDiscrAgainstMuons" , "Tau Discr. against Muons", 2, -0.5, +1.5);
00599   //hTauJetCharge_ = dqmStore_->book1D("TauJetCharge" , "Q_{#tau-Jet}", 11, -5.5, +5.5);
00600   hTauJetNumSignalTracks_ = dqmStore_->book1D("TauJetNumSignalTracks" , "Num. Tau signal Cone Tracks", 20, -0.5, +19.5);
00601   hTauJetNumIsoTracks_ = dqmStore_->book1D("TauJetNumIsoTracks" , "Num. Tau isolation Cone Tracks", 20, -0.5, +19.5);
00602   
00603   hVisMass_ = dqmStore_->book1D("VisMass", "#mu + #tau-Jet visible Mass", 20, 0., 120.);
00604 hVisMassFinal_ = dqmStore_->book1D("VisMassFinal", "#mu + #tau-Jet visible final Mass", 20, 0., 120.);
00605   //hMtMuCaloMEt_ = dqmStore_->book1D("MtMuCaloMEt", "#mu + E_{T}^{miss} (Calo) transverse Mass", 20, 20., 120.);
00606   hMtMuPFMEt_ = dqmStore_->book1D("MtMuPFMEt", "#mu + E_{T}^{miss} (PF) transverse Mass", 20, 0., 120.);
00607   //hPzetaCaloMEt_ = dqmStore_->book1D("PzetaCaloMEt", "P_{#zeta} - 1.5*P_{#zeta}^{vis} (Calo)", 20, -40., 40.);
00608   //hPzetaPFMEt_ = dqmStore_->book1D("PzetaPFMEt", "P_{#zeta} - 1.5*P_{#zeta}^{vis} (PF)", 20, -40., 40.);
00609   hMuTauAcoplanarity_ = dqmStore_->book1D("MuTauAcoplanarity", "#Delta #phi_{#mu #tau-Jet}", 20, -TMath::Pi(), +TMath::Pi());
00610   //hMuTauCharge_ = dqmStore_->book1D("MuTauCharge" , "Q_{#mu + #tau-Jet}", 11, -5.5, +5.5);
00611 hMuTauDeltaR_ = dqmStore_->book1D("MuTauDeltaR", "#Delta R_{#mu #tau-Jet}", 20, 0, 5);
00612   //hVertexChi2_ = dqmStore_->book1D("VertexChi2", "Event Vertex #chi^{2} / n.d.o.f.", 20, 0., 2.0);
00613   hVertexZ_ = dqmStore_->book1D("VertexZ", "Event Vertex z-Position", 20, -25., +25.);
00614   //hVertexD0_ = dqmStore_->book1D("VertexD0", "Event Vertex d_{0}", 20, -0.0001, 0.05);
00615 
00616   hCaloMEtPt_ = dqmStore_->book1D("CaloMEtPt", "E_{T}^{miss} (Calo)", 20, 0., 100.);
00617   //hCaloMEtPhi_ = dqmStore_->book1D("CaloMEtPhi", "#phi^{miss} (Calo)", 20, -TMath::Pi(), +TMath::Pi());
00618 
00619   hPFMEtPt_ = dqmStore_->book1D("PFMEtPt", "E_{T}^{miss} (PF)", 20, 0., 100.);
00620   //hPFMEtPhi_ = dqmStore_->book1D("PFMEtPhi", "#phi^{miss} (PF)", 20, -TMath::Pi(), +TMath::Pi());
00621 
00622   hCutFlowSummary_ = dqmStore_->book1D("CutFlowSummary", "Cut-flow Summary", 11, 0.5, 11.5);
00623   hCutFlowSummary_->setBinLabel(kPassedPreselection, "Preselection");
00624   hCutFlowSummary_->setBinLabel(kPassedTrigger, "HLT");
00625   hCutFlowSummary_->setBinLabel(kPassedMuonId, "#mu ID");
00626   hCutFlowSummary_->setBinLabel(kPassedMuonTrackIso, "#mu Trk Iso.");
00627   hCutFlowSummary_->setBinLabel(kPassedMuonEcalIso, "#mu Ecal Iso.");
00628   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrack, "#tau lead. Track");
00629   hCutFlowSummary_->setBinLabel(kPassedTauLeadTrackPt, "#tau lead. Track P_{T}");
00630   hCutFlowSummary_->setBinLabel(kPassedTauTrackIso, "#tau Track Iso.");
00631   hCutFlowSummary_->setBinLabel(kPassedTauEcalIso, "#tau Ecal Iso.");
00632   hCutFlowSummary_->setBinLabel(kPassedTauDiscrAgainstMuons, "#tau anti-#mu Discr.");
00633   //hCutFlowSummary_->setBinLabel(kPassedMuonCombIso, "#mu Comb Iso.");
00634 
00635 hCutFlowSummary_->setBinLabel(kPassedDeltaR, "#DeltaR(#mu,#tau) ");
00636 }
00637 
00638 void EwkMuTauHistManager::fillHistograms(const edm::Event& evt, const edm::EventSetup& es)
00639 {
00640   if ( cfgError_ ) return;
00641 
00642   //-----------------------------------------------------------------------------
00643   // access event-level information
00644   //-----------------------------------------------------------------------------
00645 
00646   bool readError = false;
00647 
00648 //--- get decision of high-level trigger for the event
00649   edm::Handle<edm::TriggerResults> hltDecision;
00650   readEventData(evt, triggerResultsSource_, hltDecision, numWarningsTriggerResults_, maxNumWarnings_, 
00651                 readError, "Failed to access Trigger results");
00652   if ( readError ) return;
00653   
00654   const edm::TriggerNames & triggerNames = evt.triggerNames(*hltDecision);
00655    
00656   bool isTriggered = false;
00657   for ( vstring::const_iterator hltPath = hltPaths_.begin();
00658         hltPath != hltPaths_.end(); ++hltPath ) {
00659     unsigned int index = triggerNames.triggerIndex(*hltPath);
00660     if ( index < triggerNames.size() ) {
00661       if ( hltDecision->accept(index) ) isTriggered = true;
00662     } else {
00663       if ( numWarningsHLTpath_ < maxNumWarnings_ || maxNumWarnings_ == -1 ) 
00664         edm::LogWarning ("EwkMuTauHistManager") << " Undefined HLT path = " << (*hltPath) << " !!";
00665       ++numWarningsHLTpath_;
00666       continue;
00667     }
00668   }
00669   
00670 //--- get reconstructed primary event vertex of the event
00671 //   (take as "the" primary event vertex the first entry in the collection
00672 //    of vertex objects, corresponding to the vertex associated to the highest Pt sum of tracks)
00673   edm::Handle<reco::VertexCollection> vertexCollection;
00674   readEventData(evt, vertexSource_, vertexCollection, numWarningsVertex_, maxNumWarnings_,
00675                 readError, "Failed to access Vertex collection");
00676   if ( readError ) return;
00677 
00678   const reco::Vertex* theEventVertex = ( vertexCollection->size() > 0 ) ? &(vertexCollection->at(0)) : 0;
00679 
00680 //--- get beam-spot (expected vertex position) for the event
00681   edm::Handle<reco::BeamSpot> beamSpot;
00682   readEventData(evt, beamSpotSource_, beamSpot, numWarningsBeamSpot_, maxNumWarnings_,
00683                 readError, "Failed to access Beam-spot");
00684   if ( readError ) return;
00685   
00686 //--- get collections of reconstructed muons from the event
00687   edm::Handle<reco::MuonCollection> muons;
00688   readEventData(evt, muonSource_, muons, numWarningsMuon_, maxNumWarnings_,
00689                 readError, "Failed to access Muon collection");
00690   if ( readError ) return;
00691 
00692   const reco::Muon* theMuon = getTheMuon(*muons, muonEtaCut_, muonPtCut_);
00693 
00694   double theMuonTrackIsoPt = 1.e+3;
00695   double theMuonEcalIsoPt = 1.e+3;
00696 double theMuonCombIsoPt = 1.e+3;
00697 
00698   if ( theMuon ) {
00699     theMuonTrackIsoPt = theMuon->isolationR05().sumPt;
00700     //mu.isolationR05().emEt + mu.isolationR05().hadEt + mu.isolationR05().sumPt
00701     theMuonEcalIsoPt = theMuon->isolationR05().emEt;
00702 
00703     if ( muonIsoMode_ == kRelativeIso && theMuon->pt() > 0. ) {
00704       theMuonTrackIsoPt /= theMuon->pt();
00705       theMuonEcalIsoPt /= theMuon->pt();
00706       theMuonCombIsoPt=(theMuon->isolationR05().sumPt+theMuon->isolationR05().emEt)/theMuon->pt();
00707       //std::cout<<"Rel Iso ="<<theMuonCombIsoPt<<std::endl;
00708     }
00709   }
00710 
00711 //--- get collections of reconstructed tau-jets from the event
00712   edm::Handle<reco::PFTauCollection> tauJets;
00713   readEventData(evt, tauJetSource_, tauJets, numWarningsTauJet_, maxNumWarnings_,
00714                 readError, "Failed to access Tau-jet collection");
00715   if ( readError ) return;
00716 
00717 //--- get collections of tau-jet discriminators for those tau-jets
00718   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackFinding;
00719   readEventData(evt, tauDiscrByLeadTrackFinding_, tauDiscrByLeadTrackFinding, numWarningsTauDiscrByLeadTrackFinding_, maxNumWarnings_,
00720                 readError, "Failed to access collection of pf. Tau discriminators by leading Track finding");
00721   edm::Handle<reco::PFTauDiscriminator> tauDiscrByLeadTrackPtCut;
00722   readEventData(evt, tauDiscrByLeadTrackPtCut_, tauDiscrByLeadTrackPtCut, numWarningsTauDiscrByLeadTrackPtCut_, maxNumWarnings_,
00723                 readError, "Failed to access collection of pf. Tau discriminators by leading Track Pt cut");
00724   edm::Handle<reco::PFTauDiscriminator> tauDiscrByTrackIso;
00725   readEventData(evt, tauDiscrByTrackIso_, tauDiscrByTrackIso, numWarningsTauDiscrByTrackIso_, maxNumWarnings_,
00726                 readError, "Failed to access collection of pf. Tau discriminators by Track isolation");
00727   edm::Handle<reco::PFTauDiscriminator> tauDiscrByEcalIso;
00728   readEventData(evt, tauDiscrByTrackIso_, tauDiscrByEcalIso, numWarningsTauDiscrByEcalIso_, maxNumWarnings_,
00729                 readError, "Failed to access collection of pf. Tau discriminators by ECAL isolation");
00730   edm::Handle<reco::PFTauDiscriminator> tauDiscrAgainstMuons;
00731   readEventData(evt, tauDiscrAgainstMuons_, tauDiscrAgainstMuons, numWarningsTauDiscrAgainstMuons_, maxNumWarnings_,
00732                 readError, "Failed to access collection of pf. Tau discriminators against Muons");
00733   if ( readError ) return;
00734 
00735   int theTauJetIndex = -1;
00736   const reco::PFTau* theTauJet = getTheTauJet(*tauJets, tauJetEtaCut_, tauJetPtCut_, theTauJetIndex);
00737 
00738   double theTauDiscrByLeadTrackFinding = -1.;
00739   double theTauDiscrByLeadTrackPtCut = -1.;
00740   double theTauDiscrByTrackIso = -1.;
00741   double theTauDiscrByEcalIso = -1.;
00742   double theTauDiscrAgainstMuons = -1.;
00743   if ( theTauJetIndex != -1 ) {
00744     reco::PFTauRef theTauJetRef(tauJets, theTauJetIndex);
00745     theTauDiscrByLeadTrackFinding = (*tauDiscrByLeadTrackFinding)[theTauJetRef];
00746     theTauDiscrByLeadTrackPtCut = (*tauDiscrByLeadTrackPtCut)[theTauJetRef];
00747     theTauDiscrByTrackIso = (*tauDiscrByTrackIso)[theTauJetRef];
00748     theTauDiscrByEcalIso = (*tauDiscrByEcalIso)[theTauJetRef];
00749     theTauDiscrAgainstMuons = (*tauDiscrAgainstMuons)[theTauJetRef];
00750   }
00751 
00752 //--- get missing transverse momentum
00753 //    measured by calorimeters/reconstructed by particle-flow algorithm
00754   edm::Handle<reco::CaloMETCollection> caloMEtCollection;
00755   readEventData(evt, caloMEtSource_, caloMEtCollection, numWarningsCaloMEt_, maxNumWarnings_,
00756                 readError, "Failed to access calo. MET collection");
00757   if ( readError ) return;
00758 
00759   const reco::CaloMET& caloMEt = caloMEtCollection->at(0);
00760   
00761   edm::Handle<reco::PFMETCollection> pfMEtCollection;
00762   readEventData(evt, pfMEtSource_, pfMEtCollection, numWarningsPFMEt_, maxNumWarnings_,
00763                 readError, "Failed to access pf. MET collection");
00764   if ( readError ) return;
00765 
00766   const reco::PFMET& pfMEt = pfMEtCollection->at(0);
00767 
00768   if ( !(theMuon && theTauJet && theTauJetIndex != -1) ) return;
00769 
00770   //-----------------------------------------------------------------------------
00771   // compute EWK tau analysis specific quantities
00772   //-----------------------------------------------------------------------------
00773 
00774   double dPhiMuTau = calcDeltaPhi(theMuon->phi(), theTauJet->phi());
00775   //double dRMuTau = calcDeltaR(theMuon->p4(), theTauJet->p4());
00776   double dRMuTau = fabs(ROOT::Math::VectorUtil::DeltaR(theMuon->p4(),theTauJet->p4()));
00777   double mMuTau = (theMuon->p4() + theTauJet->p4()).M();
00778 
00779   //double mtMuCaloMEt = calcMt(theMuon->px(), theMuon->px(), caloMEt.px(), caloMEt.py());
00780   double mtMuPFMEt = calcMt(theMuon->px(), theMuon->px(), pfMEt.px(), pfMEt.py());
00781 
00782   //double pZetaCaloMEt = calcPzeta(theMuon->p4(), theTauJet->p4(), caloMEt.px(), caloMEt.py());
00783   //double pZetaPFMEt = calcPzeta(theMuon->p4(), theTauJet->p4(), pfMEt.px(), pfMEt.py());
00784 
00785   //-----------------------------------------------------------------------------
00786   // apply selection criteria; fill histograms
00787   //-----------------------------------------------------------------------------
00788 
00789 //--- fill muon multiplicity histogram
00790   unsigned numGlobalMuons = 0;
00791   for ( reco::MuonCollection::const_iterator muon = muons->begin();
00792         muon != muons->end(); ++muon ) {
00793     if ( muon->isGlobalMuon() ) {
00794       ++numGlobalMuons;
00795     }
00796   }
00797 
00798   //hNumGlobalMuons_->Fill(numGlobalMuons);
00799 
00800   ++numEventsAnalyzed_;
00801 
00802   bool isSelected = false;
00803   int cutFlowStatus = -1;
00804 
00805   //if ( muonIsoMode_ == kAbsoluteIso){
00806   if ( mMuTau > visMassCut_ ) {
00807     cutFlowStatus = kPassedPreselection;
00808   }
00809   if ( cutFlowStatus == kPassedPreselection && (isTriggered || hltPaths_.size() == 0) ) {
00810     cutFlowStatus = kPassedTrigger;
00811   }
00812   if ( cutFlowStatus == kPassedTrigger && (theMuon->isGlobalMuon()||theMuon->isTrackerMuon()) ) {
00813     cutFlowStatus = kPassedMuonId;
00814    
00815   }
00816  
00817   if ( cutFlowStatus == kPassedMuonId && (theTauDiscrByLeadTrackFinding > 0.5) && (theTauJet->eta() < tauJetEtaCut_) && (theTauJet->pt() > tauJetPtCut_)  ) {
00818     cutFlowStatus = kPassedTauLeadTrack;
00819     
00820   }
00821   if ( cutFlowStatus == kPassedTauLeadTrack && theTauDiscrByLeadTrackPtCut > 0.5 ) {
00822     cutFlowStatus = kPassedTauLeadTrackPt;
00823     //hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());
00824   }
00825   if ( cutFlowStatus == kPassedTauLeadTrackPt && theTauDiscrAgainstMuons  > 0.5 ) {
00826     cutFlowStatus = kPassedTauDiscrAgainstMuons;
00827     //hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());
00828   }
00829   if ( cutFlowStatus == kPassedTauDiscrAgainstMuons && dRMuTau > deltaRCut_) {
00830     cutFlowStatus = kPassedDeltaR;
00831     //hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);
00832 
00833 
00834 
00835 hMuonPt_->Fill(theMuon->pt());
00836     hMuonEta_->Fill(theMuon->eta());
00837     hMuonPhi_->Fill(theMuon->phi());
00838 
00839     hTauJetPt_->Fill(theTauJet->pt());
00840     hTauJetEta_->Fill(theTauJet->eta());
00841     hTauJetPhi_->Fill(theTauJet->phi());
00842 
00843     //hTauJetCharge_->Fill(theTauJet->charge());
00844     if ( theTauJet->signalTracks().isAvailable()    ) hTauJetNumSignalTracks_->Fill(theTauJet->signalTracks().size());
00845     if ( theTauJet->isolationTracks().isAvailable() ) hTauJetNumIsoTracks_->Fill(theTauJet->isolationTracks().size());
00846   
00847     hVisMass_->Fill(mMuTau);
00848     //hMtMuCaloMEt_->Fill(mtMuCaloMEt);
00849     hMtMuPFMEt_->Fill(mtMuPFMEt);
00850     //hPzetaCaloMEt_->Fill(pZetaCaloMEt);
00851     //hPzetaPFMEt_->Fill(pZetaPFMEt);
00852     hMuTauAcoplanarity_->Fill(dPhiMuTau);
00853 hMuTauDeltaR_->Fill(dRMuTau);
00854     //hMuTauCharge_->Fill(theMuon->charge() + theTauJet->charge());
00855 
00856     if ( theEventVertex ) {
00857       //hVertexChi2_->Fill(theEventVertex->normalizedChi2());
00858       hVertexZ_->Fill(theEventVertex->z());
00859       //hVertexD0_->Fill(getVertexD0(*theEventVertex, *beamSpot));
00860     }
00861     
00862     hCaloMEtPt_->Fill(caloMEt.pt());
00863     //hCaloMEtPhi_->Fill(caloMEt.phi());
00864 
00865     hPFMEtPt_->Fill(pfMEt.pt());
00866     //hPFMEtPhi_->Fill(pfMEt.phi());
00867     hMuonTrackIsoPt_->Fill(theMuonTrackIsoPt);
00868 hMuonEcalIsoPt_->Fill(theMuonEcalIsoPt);
00869 hMuonCombIsoPt_->Fill(theMuonCombIsoPt);
00870 //hMuonCombIsoPt_->Fill((theMuonTrackIsoPt+theMuonEcalIsoPt)/theMuon->pt());
00871 
00872 // std::cout<<"Rel Iso Hist = "<<(theMuonTrackIsoPt+theMuonEcalIsoPt)/theMuon->pt()<<std::endl;
00873 hTauEcalIsoPt_->Fill(theTauJet->isolationPFGammaCandsEtSum());
00874 hTauTrackIsoPt_->Fill(theTauJet->isolationPFChargedHadrCandsPtSum());
00875 hTauDiscrAgainstMuons_->Fill(theTauDiscrAgainstMuons);
00876 if ( theTauJet->leadTrack().isAvailable() ) hTauLeadTrackPt_->Fill(theTauJet->leadTrack()->pt());
00877 
00878   }
00879 
00880   
00881 
00882   if ( (cutFlowStatus == kPassedDeltaR) && (((theMuonTrackIsoPt < muonTrackIsoCut_)&&(muonIsoMode_ == kAbsoluteIso))||((1>0)&&(muonIsoMode_ == kRelativeIso))) ) {
00883     cutFlowStatus = kPassedMuonTrackIso;
00884     //isSelected = true;
00885   }
00886   if ( cutFlowStatus == kPassedMuonTrackIso && (((theMuonEcalIsoPt < muonEcalIsoCut_)&&(muonIsoMode_ == kAbsoluteIso))||((theMuonCombIsoPt < muonCombIsoCut_)&&(muonIsoMode_ == kRelativeIso))) ) {
00887   cutFlowStatus = kPassedMuonEcalIso;
00888     //isSelected = true;
00889   }
00890 
00891 if ( cutFlowStatus == kPassedMuonEcalIso && theTauDiscrByTrackIso > 0.5 )
00892   {cutFlowStatus = kPassedTauTrackIso;
00893 
00894   }
00895 
00896 if ( cutFlowStatus == kPassedTauTrackIso && theTauDiscrByEcalIso > 0.5 ) {
00897 cutFlowStatus = kPassedTauEcalIso;
00898 isSelected = true;
00899 
00900  }
00901 
00902   for ( int iCut = 1; iCut <= cutFlowStatus; ++iCut ) {
00903     hCutFlowSummary_->Fill(iCut);
00904   }
00905 
00906   
00907 
00908 
00909 
00910 
00911 
00912 
00913 for ( int iCut = 1; iCut <= cutFlowStatus; ++iCut ) {
00914      hCutFlowSummary_->Fill(iCut);
00915    }
00916 
00917 
00918 
00919 //     }
00920     
00921 
00922 
00923 
00924   if ( isSelected ) 
00925     {
00926 hVisMassFinal_->Fill(mMuTau);
00927 ++numEventsSelected_;
00928     }
00929 }
00930 
00931 void EwkMuTauHistManager::finalizeHistograms()
00932 {
00933   edm::LogInfo ("EwkMuTauHistManager") 
00934     << "Filter-Statistics Summary:" << std::endl
00935     << " Events analyzed = " << numEventsAnalyzed_ << std::endl
00936     << " Events selected = " << numEventsSelected_;
00937   if ( numEventsAnalyzed_ > 0 ) {
00938     double eff = numEventsSelected_/(double)numEventsAnalyzed_;
00939     edm::LogInfo ("") 
00940       << "Overall efficiency = " << std::setprecision(4) << eff*100. 
00941       << " +/- " << std::setprecision(4) << TMath::Sqrt(eff*(1 - eff)/numEventsAnalyzed_)*100. << ")%";
00942   }
00943 }
00944 
00945 //-------------------------------------------------------------------------------
00946 // common auxiliary functions used by different channels
00947 //-------------------------------------------------------------------------------
00948 
00949 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00950 
00951 #include <TMath.h>
00952 
00953 int getIsoMode(const std::string& isoMode_string, int& error)
00954 {
00955   int isoMode_int;
00956   if ( isoMode_string == "absoluteIso" ) {
00957     isoMode_int = kAbsoluteIso;
00958   } else if ( isoMode_string == "relativeIso" ) {
00959     isoMode_int = kRelativeIso;
00960   } else { 
00961     edm::LogError ("getIsoMode") << " Failed to decode isoMode string = " << isoMode_string << " !!";
00962     isoMode_int = kUndefinedIso;
00963     error = 1;
00964   }
00965   return isoMode_int;
00966 }
00967 
00968 //
00969 //-----------------------------------------------------------------------------------------------------------------------
00970 //
00971 
00972 double calcDeltaPhi(double phi1, double phi2) 
00973 {
00974   double deltaPhi = phi1 - phi2;
00975 
00976   if ( deltaPhi < 0. ) deltaPhi = -deltaPhi;
00977 
00978   if ( deltaPhi > TMath::Pi() ) deltaPhi = 2*TMath::Pi() - deltaPhi;
00979 
00980   return deltaPhi;
00981 }
00982 
00983 double calcMt(double px1, double py1, double px2, double py2) 
00984 {
00985   double pt1 = TMath::Sqrt(px1*px1 + py1*py1);
00986   double pt2 = TMath::Sqrt(px2*px2 + py2*py2);
00987 
00988   double p1Dotp2 = px1*px2 + py1*py2;
00989   double cosAlpha = p1Dotp2/(pt1*pt2);
00990 
00991   return TMath::Sqrt(2*pt1*pt2*(1 - cosAlpha));
00992 }
00993 
00994 double calcPzeta(const reco::Candidate::LorentzVector& p1,const reco::Candidate::LorentzVector& p2, double pxMEt, double pyMEt)
00995 {
00996   double cosPhi1 = cos(p1.phi());
00997   double sinPhi1 = sin(p1.phi());
00998   double cosPhi2 = cos(p2.phi());
00999   double sinPhi2 = sin(p2.phi());
01000   double zetaX = cosPhi1 + cosPhi2;
01001   double zetaY = sinPhi1 + sinPhi2;
01002   double zetaR = TMath::Sqrt(zetaX*zetaX + zetaY*zetaY);
01003   if ( zetaR > 0. ) {
01004     zetaX /= zetaR;
01005     zetaY /= zetaR;
01006   }
01007 
01008   double pxVis = p1.px() + p2.px();
01009   double pyVis = p1.py() + p2.py();
01010   double pZetaVis = pxVis*zetaX + pyVis*zetaY;
01011 
01012   double px = pxVis + pxMEt;
01013   double py = pyVis + pyMEt;
01014   double pZeta = px*zetaX + py*zetaY;
01015 
01016   return pZeta - 1.5*pZetaVis;
01017 }
01018 
01019 //
01020 //-----------------------------------------------------------------------------------------------------------------------
01021 //
01022 
01023 bool passesElectronPreId(const reco::GsfElectron& electron)
01024 {
01025   if ( (TMath::Abs(electron.eta()) < 1.479 || TMath::Abs(electron.eta()) > 1.566) && // cut ECAL barrel/endcap crack
01026        electron.deltaPhiSuperClusterTrackAtVtx() < 0.8 &&
01027        electron.deltaEtaSuperClusterTrackAtVtx() < 0.01 &&
01028        electron.sigmaIetaIeta() < 0.03 ) {
01029     return true;
01030   } else {
01031     return false;
01032   }
01033 }
01034 
01035 bool passesElectronId(const reco::GsfElectron& electron)
01036 {
01037   if ( passesElectronPreId(electron) &&
01038        ((TMath::Abs(electron.eta()) > 1.566 && // electron reconstructed in ECAL endcap
01039          electron.sigmaEtaEta() < 0.03 && electron.hcalOverEcal() < 0.05 && 
01040          TMath::Abs(electron.deltaEtaSuperClusterTrackAtVtx()) < 0.009 && 
01041          TMath::Abs(electron.deltaPhiSuperClusterTrackAtVtx()) < 0.7 ) ||
01042         (TMath::Abs(electron.eta()) < 1.479 && // electron reconstructed in ECAL barrel
01043          electron.sigmaEtaEta() < 0.01 && electron.hcalOverEcal() < 0.12 && 
01044          TMath::Abs(electron.deltaEtaSuperClusterTrackAtVtx()) < 0.007 && 
01045          TMath::Abs(electron.deltaPhiSuperClusterTrackAtVtx()) < 0.8)) ) {
01046     return true;
01047   } else {
01048     return false;
01049   }
01050 }
01051 
01052 //
01053 //-----------------------------------------------------------------------------------------------------------------------
01054 //
01055 
01056 const reco::GsfElectron* getTheElectron(const reco::GsfElectronCollection& electrons, double electronEtaCut, double electronPtCut)
01057 {
01058   const reco::GsfElectron* theElectron = 0;
01059   
01060   for ( reco::GsfElectronCollection::const_iterator electron = electrons.begin();
01061         electron != electrons.end(); ++electron ) {
01062     if ( TMath::Abs(electron->eta()) < electronEtaCut && electron->pt() > electronPtCut && passesElectronPreId(*electron) ) {
01063       if ( theElectron == 0 || electron->pt() > theElectron->pt() ) theElectron = &(*electron);
01064     }
01065   }
01066   
01067   return theElectron;
01068 }
01069 
01070 
01071 
01072 
01073 const reco::Muon* getTheMuon(const reco::MuonCollection& muons, double muonEtaCut, double muonPtCut)
01074 {
01075   const reco::Muon* theMuon = 0;
01076 
01077   for ( reco::MuonCollection::const_iterator muon = muons.begin();
01078         muon != muons.end(); ++muon ) {
01079     if ( TMath::Abs(muon->eta()) < muonEtaCut && muon->pt() > muonPtCut ) {
01080       if ( theMuon == 0 || muon->pt() > theMuon->pt() ) theMuon = &(*muon);
01081     }
01082   }
01083 
01084   return theMuon;
01085 }
01086 
01087 const reco::PFTau* getTheTauJet(const reco::PFTauCollection& tauJets, double tauJetEtaCut, double tauJetPtCut, int& theTauJetIndex)
01088 {
01089   const reco::PFTau* theTauJet = 0;
01090   theTauJetIndex = -1;
01091 
01092   int numTauJets = tauJets.size();
01093   for ( int iTauJet = 0; iTauJet < numTauJets; ++iTauJet ) {
01094     const reco::PFTau& tauJet = tauJets.at(iTauJet);
01095 
01096     if ( fabs(tauJet.eta()) < tauJetEtaCut && tauJet.pt() > tauJetPtCut ) {
01097       if ( theTauJet == 0 || tauJet.pt() > theTauJet->pt() ) {
01098         theTauJet = &tauJet;
01099         theTauJetIndex = iTauJet;
01100       }
01101     }
01102   }
01103 
01104   return theTauJet;
01105 }
01106 
01107 //
01108 //-----------------------------------------------------------------------------------------------------------------------
01109 //
01110 
01111 double getVertexD0(const reco::Vertex& vertex, const reco::BeamSpot& beamSpot)
01112 {
01113   double dX = vertex.x() - beamSpot.x0();
01114   double dY = vertex.y() - beamSpot.y0();
01115   return TMath::Sqrt(dX*dX + dY*dY);
01116 }
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/DQM/Physics/src/EwkTauDQM.cc
