Reference Manual

Go to the documentation of this file.
00001 /* This Class Header */
00002 #include "DQMOffline/Muon/src/EfficiencyAnalyzer.h"
00003 
00004 /* Collaborating Class Header */
00005 #include "FWCore/Framework/interface/MakerMacros.h"
00006 #include "FWCore/Framework/interface/Frameworkfwd.h"
00007 #include "FWCore/Framework/interface/Event.h"
00008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00009 #include "FWCore/Framework/interface/ESHandle.h"
00010 #include "FWCore/Utilities/interface/Exception.h"
00011 #include "DataFormats/MuonReco/interface/Muon.h"
00012 #include "DataFormats/MuonReco/interface/MuonFwd.h" 
00013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00014 #include "DataFormats/Math/interface/deltaR.h"
00015 #include "DataFormats/VertexReco/interface/Vertex.h"
00016 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00017 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00018 
00019 
00020 
00021 using namespace edm;
00022 
00023 #include "TLorentzVector.h"
00024 #include "TFile.h"
00025 #include <vector>
00026 #include "math.h"
00027 #include <algorithm>
00028 
00029 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00030 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00031 
00032 
00033 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00034 #include "DataFormats/TrackReco/interface/Track.h"
00035 
00036 /* C++ Headers */
00037 #include <iostream>
00038 #include <fstream>
00039 #include <cmath>
00040 using namespace std;
00041 using namespace edm;
00042 
00043 EfficiencyAnalyzer::EfficiencyAnalyzer(const edm::ParameterSet& pSet, MuonServiceProxy *theService):MuonAnalyzerBase(theService){ 
00044   parameters = pSet;
00045 }
00046 
00047 EfficiencyAnalyzer::~EfficiencyAnalyzer() { }
00048 
00049 void EfficiencyAnalyzer::beginJob(DQMStore * dbe) {
00050 #ifdef DEBUG
00051   cout << "[EfficiencyAnalyzer] Parameters initialization" <<endl;
00052 #endif
00053   metname = "EfficiencyAnalyzer";
00054   LogTrace(metname)<<"[EfficiencyAnalyzer] Parameters initialization";
00055   dbe->setCurrentFolder("Muons/EfficiencyAnalyzer");  
00056   
00057   theMuonCollectionLabel = parameters.getParameter<edm::InputTag>("MuonCollection");
00058   theTrackCollectionLabel = parameters.getParameter<edm::InputTag>("TrackCollection");
00059 
00060 
00061   //Vertex requirements
00062   _doPVCheck = parameters.getParameter<bool>("doPrimaryVertexCheck");
00063   vertexTag  = parameters.getParameter<edm::InputTag>("vertexLabel");
00064   bsTag  = parameters.getParameter<edm::InputTag>("bsLabel");
00065 
00066   ptBin_ = parameters.getParameter<int>("ptBin");
00067   ptMin_ = parameters.getParameter<double>("ptMin");
00068   ptMax_ = parameters.getParameter<double>("ptMax");
00069 
00070   etaBin_ = parameters.getParameter<int>("etaBin");
00071   etaMin_ = parameters.getParameter<double>("etaMin");
00072   etaMax_ = parameters.getParameter<double>("etaMax");
00073 
00074   phiBin_ = parameters.getParameter<int>("phiBin");
00075   phiMin_ = parameters.getParameter<double>("phiMin");
00076   phiMax_ = parameters.getParameter<double>("phiMax");
00077 
00078   vtxBin_ = parameters.getParameter<int>("vtxBin");
00079   vtxMin_ = parameters.getParameter<double>("vtxMin");
00080   vtxMax_ = parameters.getParameter<double>("vtxMax");
00081 
00082   test_TightMu_Minv  = dbe->book1D("test_TightMu_Minv"  ,"Minv",50,70,110);
00083 
00084   h_allProbes_pt = dbe->book1D("allProbes_pt","All Probes Pt", ptBin_, ptMin_, ptMax_);
00085   h_allProbes_EB_pt = dbe->book1D("allProbes_EB_pt","Barrel: all Probes Pt", ptBin_, ptMin_, ptMax_);
00086   h_allProbes_EE_pt = dbe->book1D("allProbes_EE_pt","Endcap: all Probes Pt", ptBin_, ptMin_, ptMax_);
00087   h_allProbes_eta = dbe->book1D("allProbes_eta","All Probes Eta", etaBin_, etaMin_, etaMax_);
00088   h_allProbes_hp_eta = dbe->book1D("allProbes_hp_eta","High Pt all Probes Eta", etaBin_, etaMin_, etaMax_);
00089   h_allProbes_phi = dbe->book1D("allProbes_phi","All Probes Phi", phiBin_, phiMin_, phiMax_);
00090 
00091   h_allProbes_TightMu_pt = dbe->book1D("allProbes_TightMu_pt","All TightMu Probes Pt", ptBin_, ptMin_, ptMax_);
00092   h_allProbes_EB_TightMu_pt = dbe->book1D("allProbes_EB_TightMu_pt","Barrel: all TightMu Probes Pt", ptBin_, ptMin_, ptMax_);
00093   h_allProbes_EE_TightMu_pt = dbe->book1D("allProbes_EE_TightMu_pt","Endcap: all TightMu Probes Pt", ptBin_, ptMin_, ptMax_);
00094   h_allProbes_TightMu_nVtx = dbe->book1D("allProbes_TightMu_nVtx","All Probes (TightMu) nVtx", vtxBin_, vtxMin_, vtxMax_);
00095   h_allProbes_EB_TightMu_nVtx = dbe->book1D("allProbes_EB_TightMu_nVtx","Barrel: All Probes (TightMu) nVtx", vtxBin_, vtxMin_, vtxMax_);
00096   h_allProbes_EE_TightMu_nVtx = dbe->book1D("allProbes_EE_TightMu_nVtx","Endcap: All Probes (TightMu) nVtx", vtxBin_, vtxMin_, vtxMax_);
00097 
00098   h_passProbes_TightMu_pt = dbe->book1D("passProbes_TightMu_pt","TightMu Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
00099   h_passProbes_TightMu_EB_pt = dbe->book1D("passProbes_TightMu_EB_pt","Barrel: TightMu Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
00100   h_passProbes_TightMu_EE_pt = dbe->book1D("passProbes_TightMu_EE_pt","Endcap: TightMu Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
00101   h_passProbes_TightMu_eta = dbe->book1D("passProbes_TightMu_eta","TightMu Passing Probes #eta", etaBin_, etaMin_, etaMax_);
00102   h_passProbes_TightMu_hp_eta = dbe->book1D("passProbes_TightMu_hp_eta","High Pt TightMu Passing Probes #eta", etaBin_, etaMin_, etaMax_);
00103   h_passProbes_TightMu_phi = dbe->book1D("passProbes_TightMu_phi","TightMu Passing Probes #phi", phiBin_, phiMin_, phiMax_);
00104 
00105   h_passProbes_detIsoTightMu_pt = dbe->book1D("passProbes_detIsoTightMu_pt","detIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
00106   h_passProbes_EB_detIsoTightMu_pt = dbe->book1D("passProbes_EB_detIsoTightMu_pt","Barrel: detIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
00107   h_passProbes_EE_detIsoTightMu_pt = dbe->book1D("passProbes_EE_detIsoTightMu_pt","Endcap: detIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
00108 
00109   h_passProbes_pfIsoTightMu_pt = dbe->book1D("passProbes_pfIsoTightMu_pt","pfIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
00110   h_passProbes_EB_pfIsoTightMu_pt = dbe->book1D("passProbes_EB_pfIsoTightMu_pt","Barrel: pfIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
00111   h_passProbes_EE_pfIsoTightMu_pt = dbe->book1D("passProbes_EE_pfIsoTightMu_pt","Endcap: pfIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
00112 
00113   h_passProbes_detIsoTightMu_nVtx    = dbe->book1D("passProbes_detIsoTightMu_nVtx",    "detIsoTightMu Passing Probes nVtx (R03)",  vtxBin_, vtxMin_, vtxMax_);
00114   h_passProbes_pfIsoTightMu_nVtx     = dbe->book1D("passProbes_pfIsoTightMu_nVtx",    "pfIsoTightMu Passing Probes nVtx (R04)",  vtxBin_, vtxMin_, vtxMax_);
00115   h_passProbes_EB_detIsoTightMu_nVtx = dbe->book1D("passProbes_EB_detIsoTightMu_nVtx","Barrel: detIsoTightMu Passing Probes nVtx (R03)",  vtxBin_, vtxMin_, vtxMax_);
00116   h_passProbes_EE_detIsoTightMu_nVtx = dbe->book1D("passProbes_EE_detIsoTightMu_nVtx","Endcap: detIsoTightMu Passing Probes nVtx (R03)",  vtxBin_, vtxMin_, vtxMax_);
00117   h_passProbes_EB_pfIsoTightMu_nVtx  = dbe->book1D("passProbes_EB_pfIsoTightMu_nVtx", "Barrel: pfIsoTightMu Passing Probes nVtx (R04)",  vtxBin_, vtxMin_, vtxMax_);
00118   h_passProbes_EE_pfIsoTightMu_nVtx  = dbe->book1D("passProbes_EE_pfIsoTightMu_nVtx", "Endcap: pfIsoTightMu Passing Probes nVtx (R04)",  vtxBin_, vtxMin_, vtxMax_);
00119 
00120   
00121   // Apply deltaBeta PU corrections to the PF isolation eficiencies.
00122   
00123   h_passProbes_pfIsodBTightMu_pt = dbe->book1D("passProbes_pfIsodBTightMu_pt","pfIsoTightMu Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
00124   h_passProbes_EB_pfIsodBTightMu_pt = dbe->book1D("passProbes_EB_pfIsodBTightMu_pt","Barrel: pfIsoTightMu Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
00125   h_passProbes_EE_pfIsodBTightMu_pt = dbe->book1D("passProbes_EE_pfIsodBTightMu_pt","Endcap: pfIsoTightMu Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
00126   h_passProbes_pfIsodBTightMu_nVtx     = dbe->book1D("passProbes_pfIsodBTightMu_nVtx",    "pfIsoTightMu Passing Probes nVtx (R04) (deltaB PU correction)",  vtxBin_, vtxMin_, vtxMax_);
00127 h_passProbes_EB_pfIsodBTightMu_nVtx  = dbe->book1D("passProbes_EB_pfIsodBTightMu_nVtx", "Barrel: pfIsoTightMu Passing Probes nVtx (R04) (deltaB PU correction)",  vtxBin_, vtxMin_, vtxMax_);
00128   h_passProbes_EE_pfIsodBTightMu_nVtx  = dbe->book1D("passProbes_EE_pfIsodBTightMu_nVtx", "Endcap: pfIsoTightMu Passing Probes nVtx (R04) (deltaB PU correction)",  vtxBin_, vtxMin_, vtxMax_);
00129 
00130 
00131 
00132 #ifdef DEBUG
00133   cout << "[EfficiencyAnalyzer] Parameters initialization DONE" <<endl;
00134 #endif
00135 }
00136 
00137 void EfficiencyAnalyzer::analyze(const edm::Event & iEvent,const edm::EventSetup& iSetup) {
00138 
00139   LogTrace(metname)<<"[EfficiencyAnalyzer] Analyze the mu in different eta regions";
00140   
00141   edm::Handle<reco::MuonCollection> muons;
00142   iEvent.getByLabel(theMuonCollectionLabel, muons);
00143 
00144 
00145   edm::Handle<reco::TrackCollection> tracks;
00146   iEvent.getByLabel(theTrackCollectionLabel,tracks); 
00147 
00148 
00149   reco::BeamSpot beamSpot;
00150   Handle<reco::BeamSpot> beamSpotHandle;
00151   iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
00152   beamSpot = *beamSpotHandle;
00153 
00154 
00155   // ==========================================================
00156   //Vertex information
00157   
00158   edm::Handle<reco::VertexCollection> vertex;
00159   iEvent.getByLabel(vertexTag, vertex);
00160 
00161   _numPV = 0;
00162   bool bPrimaryVertex = true;
00163   if(_doPVCheck){ 
00164     bPrimaryVertex = false;
00165          
00166     if (!vertex.isValid()) {
00167       LogTrace(metname) << "[EfficiencyAnalyzer] Could not find vertex collection" << std::endl;
00168       bPrimaryVertex = false;
00169     }
00170     
00171     if ( vertex.isValid() ){
00172       const reco::VertexCollection& vertexCollection = *(vertex.product());
00173       int vertex_number     = vertexCollection.size();
00174       
00175       reco::VertexCollection::const_iterator v = vertexCollection.begin();
00176       for ( ; v != vertexCollection.end(); ++v) {
00177         double vertex_chi2    = v->normalizedChi2();
00178         double vertex_ndof    = v->ndof();
00179         bool   fakeVtx        = v->isFake();
00180         double vertex_Z       = v->z();
00181         
00182         if (  !fakeVtx
00183               && vertex_number >= 1
00184               && vertex_ndof   > 4
00185               && vertex_chi2   < 999
00186               && fabs(vertex_Z)< 24. ) {
00187           bPrimaryVertex = true;
00188           ++_numPV;
00189         }
00190       }
00191     }
00192   }
00193   // ==========================================================
00194 
00195 
00196   // ==========================================================
00197   // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
00198 
00199   reco::Vertex::Point posVtx;
00200   reco::Vertex::Error errVtx;
00201 
00202   unsigned int theIndexOfThePrimaryVertex = 999.;
00203  
00204   
00205   if ( vertex.isValid() ){
00206 
00207     for (unsigned int ind=0; ind<vertex->size(); ++ind) {
00208     
00209       if ( (*vertex)[ind].isValid() && !((*vertex)[ind].isFake()) ) {
00210       
00211         theIndexOfThePrimaryVertex = ind;
00212         break;
00213       }
00214     }
00215   }
00216 
00217   if (theIndexOfThePrimaryVertex<100) {
00218     posVtx = ((*vertex)[theIndexOfThePrimaryVertex]).position();
00219     errVtx = ((*vertex)[theIndexOfThePrimaryVertex]).error();
00220   }   else {
00221 
00222     LogInfo("EfficiencyAnalyzer") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
00223   
00224 
00225     edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00226     iEvent.getByLabel(bsTag,recoBeamSpotHandle);
00227     
00228     reco::BeamSpot bs = *recoBeamSpotHandle;
00229     
00230     posVtx = bs.position();
00231     errVtx(0,0) = bs.BeamWidthX();
00232     errVtx(1,1) = bs.BeamWidthY();
00233     errVtx(2,2) = bs.sigmaZ();
00234   
00235       }
00236   const reco::Vertex thePrimaryVertex(posVtx,errVtx);
00237 
00238   // ==========================================================
00239 
00240 
00241 
00242   
00243   
00244   if(!muons.isValid()) return;
00245 
00246 
00247   // Loop on muon collection
00248   TLorentzVector Mu1, Mu2;
00249 
00250   bool isMB = false;
00251   bool isME = false;
00252 
00253   for (reco::MuonCollection::const_iterator recoMu1 = muons->begin(); recoMu1!=muons->end(); ++recoMu1) {
00254     LogTrace(metname)<<"[EfficiencyAnalyzer] loop over first muons" << endl;
00255     
00256     //--- Define combined isolation
00257     reco::MuonIsolation Iso_muon = recoMu1->isolationR03();
00258     float combIso = (Iso_muon.emEt + Iso_muon.hadEt + Iso_muon.sumPt);  
00259     
00260     //--- Is Global Muon 
00261     if (!recoMu1->isGlobalMuon()) continue;
00262 
00263     // get the track combinig the information from both the Tracker and the Spectrometer
00264     reco::TrackRef recoCombinedGlbTrack1 = recoMu1->combinedMuon();    
00265     float muPt1 = recoCombinedGlbTrack1->pt();
00266     Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(), recoCombinedGlbTrack1->py(),recoCombinedGlbTrack1->pz(), recoCombinedGlbTrack1->p());
00267     
00268     
00269     //--- Define if it is a tight muon
00270     // Change the Tight muon definition by using the implemented method in: MuonSelectors.cc
00271 
00272     if (!muon::isTightMuon(*recoMu1, thePrimaryVertex)) continue;
00273   
00274     //-- is isolated muon
00275     if (muPt1 <= 15)          continue;
00276     if (combIso/muPt1 > 0.1 ) continue;
00277     
00278     for (reco::MuonCollection::const_iterator recoMu2 = muons->begin(); recoMu2!=muons->end(); ++recoMu2){ 
00279       LogTrace(metname)<<"[EfficiencyAnalyzer] loop over second muon" <<endl;
00280       if (recoMu2 == recoMu1) continue;
00281       
00282       if (recoMu2->eta() < 1.479 )  isMB = true;
00283       if (recoMu2->eta() >= 1.479 ) isME = true;
00284       
00285       //--> should we apply track quality cuts??? 
00286       Mu2.SetPxPyPzE(recoMu2->px(), recoMu2->py(), recoMu2->pz(), recoMu2->p());
00287       
00288       float Minv = (Mu1+Mu2).M();
00289       if (!recoMu2->isTrackerMuon())                     continue;
00290       if ( recoMu2->pt() < 5 )                           continue;
00291       if ( (recoMu1->charge())*(recoMu2->charge()) > 0 ) continue; 
00292       if ( Minv < 70 ||  Minv > 110 )                    continue;
00293       
00294       h_allProbes_pt->Fill(recoMu2->pt());
00295       h_allProbes_eta->Fill(recoMu2->eta());
00296       h_allProbes_phi->Fill(recoMu2->phi());
00297       
00298       if (isMB)               h_allProbes_EB_pt->Fill(recoMu2->pt());
00299       if (isME)               h_allProbes_EE_pt->Fill(recoMu2->pt());
00300       if(recoMu2->pt() > 20 ) h_allProbes_hp_eta->Fill(recoMu2->eta());
00301 
00302       test_TightMu_Minv->Fill(Minv);
00303       
00304 
00305 
00306       // Probes passing the tight muon criteria 
00307 
00308       if (!muon::isTightMuon(*recoMu2, thePrimaryVertex)) continue;
00309 
00310       h_passProbes_TightMu_pt->Fill(recoMu2->pt());
00311       h_passProbes_TightMu_eta->Fill(recoMu2->eta());
00312       h_passProbes_TightMu_phi->Fill(recoMu2->phi());
00313       
00314       if (isMB) h_passProbes_TightMu_EB_pt->Fill(recoMu2->pt());
00315       if (isME) h_passProbes_TightMu_EE_pt->Fill(recoMu2->pt());
00316       if( recoMu2->pt() > 20 ) h_passProbes_TightMu_hp_eta->Fill(recoMu2->eta());
00317       
00318       h_allProbes_TightMu_pt->Fill(recoMu2->pt());
00319       if (isMB) h_allProbes_EB_TightMu_pt->Fill(recoMu2->pt());
00320       if (isME) h_allProbes_EE_TightMu_pt->Fill(recoMu2->pt());
00321 
00322       //------- For PU monitoring -------//
00323       if (bPrimaryVertex)         h_allProbes_TightMu_nVtx->Fill(_numPV);
00324       if (bPrimaryVertex && isMB) h_allProbes_EB_TightMu_nVtx->Fill(_numPV);
00325       if (bPrimaryVertex && isME) h_allProbes_EE_TightMu_nVtx->Fill(_numPV);
00326 
00327       //-- Define det relative isolation
00328       float tkIso = recoMu2->isolationR03().sumPt;
00329       float emIso = recoMu2->isolationR03().emEt;
00330       float hadIso = recoMu2->isolationR03().hadEt + recoMu2->isolationR03().hoEt;
00331       float relDetIso = (tkIso + emIso + hadIso) /  (recoMu2->pt()); 
00332  
00333       if (relDetIso < 0.05 ) { 
00334         h_passProbes_detIsoTightMu_pt->Fill(recoMu2->pt());
00335         if (isMB) h_passProbes_EB_detIsoTightMu_pt->Fill(recoMu2->pt());
00336         if (isME) h_passProbes_EE_detIsoTightMu_pt->Fill(recoMu2->pt());
00337         
00338         if (bPrimaryVertex)         h_passProbes_detIsoTightMu_nVtx->Fill(_numPV);
00339         if (bPrimaryVertex && isMB) h_passProbes_EB_detIsoTightMu_nVtx->Fill(_numPV);
00340         if (bPrimaryVertex && isME) h_passProbes_EE_detIsoTightMu_nVtx->Fill(_numPV);
00341       }
00342       
00343       //-- Define PF relative isolation
00344       float chargedIso = recoMu2->pfIsolationR04().sumChargedHadronPt;
00345       float neutralIso = recoMu2->pfIsolationR04().sumNeutralHadronEt;
00346       float photonIso = recoMu2->pfIsolationR04().sumPhotonEt;
00347       float relPFIso = (chargedIso + neutralIso + photonIso) /  (recoMu2->pt()); 
00348       
00349       float pu = recoMu2->pfIsolationR04().sumPUPt; 
00350 
00351       float neutralphotonPUCorrected = std::max(0.0 , (neutralIso + photonIso - 0.5*pu ) ); 
00352 
00353       float relPFIsoPUCorrected = (chargedIso + neutralphotonPUCorrected) /  (recoMu2->pt()); 
00354 
00355 
00356 
00357       if (relPFIso < 0.12 ) { 
00358         h_passProbes_pfIsoTightMu_pt->Fill(recoMu2->pt());
00359         if (isMB) h_passProbes_EB_pfIsoTightMu_pt->Fill(recoMu2->pt());
00360         if (isME) h_passProbes_EE_pfIsoTightMu_pt->Fill(recoMu2->pt());
00361         
00362         if( bPrimaryVertex)         h_passProbes_pfIsoTightMu_nVtx->Fill(_numPV);
00363         if (bPrimaryVertex && isMB) h_passProbes_EB_pfIsoTightMu_nVtx->Fill(_numPV);
00364         if (bPrimaryVertex && isME) h_passProbes_EE_pfIsoTightMu_nVtx->Fill(_numPV);
00365       }
00366 
00367 
00368 
00369         // Apply deltaBeta PU corrections to the PF isolation eficiencies.
00370 
00371         if ( relPFIsoPUCorrected < 0.12 ) { 
00372 
00373         h_passProbes_pfIsodBTightMu_pt->Fill(recoMu2->pt());
00374         if (isMB) h_passProbes_EB_pfIsodBTightMu_pt->Fill(recoMu2->pt());
00375         if (isME) h_passProbes_EE_pfIsodBTightMu_pt->Fill(recoMu2->pt());
00376         
00377         if( bPrimaryVertex)         h_passProbes_pfIsodBTightMu_nVtx->Fill(_numPV);
00378         if (bPrimaryVertex && isMB) h_passProbes_EB_pfIsodBTightMu_nVtx->Fill(_numPV);
00379         if (bPrimaryVertex && isME) h_passProbes_EE_pfIsodBTightMu_nVtx->Fill(_numPV);
00380            }
00381 
00382     }
00383   }
00384 }
00385 
00386
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/DQMOffline/Muon/src/EfficiencyAnalyzer.cc
