CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EfficiencyAnalyzer.cc
Go to the documentation of this file.
1 /* This Class Header */
3 
4 /* Collaborating Class Header */
15 
16 #include "TLorentzVector.h"
17 #include "TFile.h"
18 #include <vector>
19 #include <cmath>
20 #include <algorithm>
21 
22 /* C++ Headers */
23 #include <iostream>
24 #include <fstream>
25 #include <cmath>
26 using namespace std;
27 using namespace edm;
28 
30  parameters = pSet;
31 
32  theService = new MuonServiceProxy(parameters.getParameter<ParameterSet>("ServiceParameters"));
33  theDbe = edm::Service<DQMStore>().operator->();
34 
35  // DATA
36  theMuonCollectionLabel_ = consumes<reco::MuonCollection> (parameters.getParameter<edm::InputTag>("MuonCollection"));
37  theTrackCollectionLabel_ = consumes<reco::TrackCollection> (parameters.getParameter<edm::InputTag>("TrackCollection"));
38  theVertexLabel_ = consumes<reco::VertexCollection>(parameters.getParameter<edm::InputTag>("VertexLabel"));
39  theBeamSpotLabel_ = mayConsume<reco::BeamSpot> (parameters.getParameter<edm::InputTag>("BeamSpotLabel"));
40 
41 }
42 
44  delete theService;
45 }
46 
48  metname = "EfficiencyAnalyzer";
49  LogTrace(metname)<<"[EfficiencyAnalyzer] Parameters initialization";
50 }
51 
52 void EfficiencyAnalyzer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup){
53  metname = "EfficiencyAnalyzer";
54 
55  theDbe->cd();
56  theDbe->setCurrentFolder("Muons/EfficiencyAnalyzer");
57 
58  //Vertex requirements
59  doPVCheck_ = parameters.getParameter<bool>("doPrimaryVertexCheck");
60 
61  ptBin_ = parameters.getParameter<int>("ptBin");
62  ptMin_ = parameters.getParameter<double>("ptMin");
63  ptMax_ = parameters.getParameter<double>("ptMax");
64 
65  etaBin_ = parameters.getParameter<int>("etaBin");
66  etaMin_ = parameters.getParameter<double>("etaMin");
67  etaMax_ = parameters.getParameter<double>("etaMax");
68 
69  phiBin_ = parameters.getParameter<int>("phiBin");
70  phiMin_ = parameters.getParameter<double>("phiMin");
71  phiMax_ = parameters.getParameter<double>("phiMax");
72 
73  vtxBin_ = parameters.getParameter<int>("vtxBin");
74  vtxMin_ = parameters.getParameter<double>("vtxMin");
75  vtxMax_ = parameters.getParameter<double>("vtxMax");
76 
77  test_TightMu_Minv = theDbe->book1D("test_TightMu_Minv" ,"Minv",50,70,110);
78 
79  h_allProbes_pt = theDbe->book1D("allProbes_pt","All Probes Pt", ptBin_, ptMin_, ptMax_);
80  h_allProbes_EB_pt = theDbe->book1D("allProbes_EB_pt","Barrel: all Probes Pt", ptBin_, ptMin_, ptMax_);
81  h_allProbes_EE_pt = theDbe->book1D("allProbes_EE_pt","Endcap: all Probes Pt", ptBin_, ptMin_, ptMax_);
82  h_allProbes_eta = theDbe->book1D("allProbes_eta","All Probes Eta", etaBin_, etaMin_, etaMax_);
83  h_allProbes_hp_eta = theDbe->book1D("allProbes_hp_eta","High Pt all Probes Eta", etaBin_, etaMin_, etaMax_);
84  h_allProbes_phi = theDbe->book1D("allProbes_phi","All Probes Phi", phiBin_, phiMin_, phiMax_);
85 
86  h_allProbes_TightMu_pt = theDbe->book1D("allProbes_TightMu_pt","All TightMu Probes Pt", ptBin_, ptMin_, ptMax_);
87  h_allProbes_EB_TightMu_pt = theDbe->book1D("allProbes_EB_TightMu_pt","Barrel: all TightMu Probes Pt", ptBin_, ptMin_, ptMax_);
88  h_allProbes_EE_TightMu_pt = theDbe->book1D("allProbes_EE_TightMu_pt","Endcap: all TightMu Probes Pt", ptBin_, ptMin_, ptMax_);
89  h_allProbes_TightMu_nVtx = theDbe->book1D("allProbes_TightMu_nVtx","All Probes (TightMu) nVtx", vtxBin_, vtxMin_, vtxMax_);
90  h_allProbes_EB_TightMu_nVtx = theDbe->book1D("allProbes_EB_TightMu_nVtx","Barrel: All Probes (TightMu) nVtx", vtxBin_, vtxMin_, vtxMax_);
91  h_allProbes_EE_TightMu_nVtx = theDbe->book1D("allProbes_EE_TightMu_nVtx","Endcap: All Probes (TightMu) nVtx", vtxBin_, vtxMin_, vtxMax_);
92 
93  h_passProbes_TightMu_pt = theDbe->book1D("passProbes_TightMu_pt","TightMu Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
94  h_passProbes_TightMu_EB_pt = theDbe->book1D("passProbes_TightMu_EB_pt","Barrel: TightMu Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
95  h_passProbes_TightMu_EE_pt = theDbe->book1D("passProbes_TightMu_EE_pt","Endcap: TightMu Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
96  h_passProbes_TightMu_eta = theDbe->book1D("passProbes_TightMu_eta","TightMu Passing Probes #eta", etaBin_, etaMin_, etaMax_);
97  h_passProbes_TightMu_hp_eta = theDbe->book1D("passProbes_TightMu_hp_eta","High Pt TightMu Passing Probes #eta", etaBin_, etaMin_, etaMax_);
98  h_passProbes_TightMu_phi = theDbe->book1D("passProbes_TightMu_phi","TightMu Passing Probes #phi", phiBin_, phiMin_, phiMax_);
99 
100  h_passProbes_detIsoTightMu_pt = theDbe->book1D("passProbes_detIsoTightMu_pt","detIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
101  h_passProbes_EB_detIsoTightMu_pt = theDbe->book1D("passProbes_EB_detIsoTightMu_pt","Barrel: detIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
102  h_passProbes_EE_detIsoTightMu_pt = theDbe->book1D("passProbes_EE_detIsoTightMu_pt","Endcap: detIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
103 
104  h_passProbes_pfIsoTightMu_pt = theDbe->book1D("passProbes_pfIsoTightMu_pt","pfIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
105  h_passProbes_EB_pfIsoTightMu_pt = theDbe->book1D("passProbes_EB_pfIsoTightMu_pt","Barrel: pfIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
106  h_passProbes_EE_pfIsoTightMu_pt = theDbe->book1D("passProbes_EE_pfIsoTightMu_pt","Endcap: pfIsoTightMu Passing Probes Pt", ptBin_, ptMin_, ptMax_);
107 
108  h_passProbes_detIsoTightMu_nVtx = theDbe->book1D("passProbes_detIsoTightMu_nVtx", "detIsoTightMu Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
109  h_passProbes_pfIsoTightMu_nVtx = theDbe->book1D("passProbes_pfIsoTightMu_nVtx", "pfIsoTightMu Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
110  h_passProbes_EB_detIsoTightMu_nVtx = theDbe->book1D("passProbes_EB_detIsoTightMu_nVtx","Barrel: detIsoTightMu Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
111  h_passProbes_EE_detIsoTightMu_nVtx = theDbe->book1D("passProbes_EE_detIsoTightMu_nVtx","Endcap: detIsoTightMu Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
112  h_passProbes_EB_pfIsoTightMu_nVtx = theDbe->book1D("passProbes_EB_pfIsoTightMu_nVtx", "Barrel: pfIsoTightMu Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
113  h_passProbes_EE_pfIsoTightMu_nVtx = theDbe->book1D("passProbes_EE_pfIsoTightMu_nVtx", "Endcap: pfIsoTightMu Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
114 
115 
116  // Apply deltaBeta PU corrections to the PF isolation eficiencies.
117 
118  h_passProbes_pfIsodBTightMu_pt = theDbe->book1D("passProbes_pfIsodBTightMu_pt","pfIsoTightMu Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
119  h_passProbes_EB_pfIsodBTightMu_pt = theDbe->book1D("passProbes_EB_pfIsodBTightMu_pt","Barrel: pfIsoTightMu Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
120  h_passProbes_EE_pfIsodBTightMu_pt = theDbe->book1D("passProbes_EE_pfIsodBTightMu_pt","Endcap: pfIsoTightMu Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
121  h_passProbes_pfIsodBTightMu_nVtx = theDbe->book1D("passProbes_pfIsodBTightMu_nVtx", "pfIsoTightMu Passing Probes nVtx (R04) (deltaB PU correction)", vtxBin_, vtxMin_, vtxMax_);
122 h_passProbes_EB_pfIsodBTightMu_nVtx = theDbe->book1D("passProbes_EB_pfIsodBTightMu_nVtx", "Barrel: pfIsoTightMu Passing Probes nVtx (R04) (deltaB PU correction)", vtxBin_, vtxMin_, vtxMax_);
123  h_passProbes_EE_pfIsodBTightMu_nVtx = theDbe->book1D("passProbes_EE_pfIsodBTightMu_nVtx", "Endcap: pfIsoTightMu Passing Probes nVtx (R04) (deltaB PU correction)", vtxBin_, vtxMin_, vtxMax_);
124 
125 
126 
127 #ifdef DEBUG
128  cout << "[EfficiencyAnalyzer] Parameters initialization DONE" <<endl;
129 #endif
130 }
131 
133 
134  LogTrace(metname)<<"[EfficiencyAnalyzer] Analyze the mu in different eta regions";
135  theService->update(iSetup);
136  // ==========================================================
137  // BEGIN READ DATA:
138  // Muon information
140  iEvent.getByToken(theMuonCollectionLabel_, muons);
141 
142  // Tracks information
144  iEvent.getByToken(theTrackCollectionLabel_,tracks);
145 
146  //Vertex information
148  iEvent.getByToken(theVertexLabel_, vertex);
149  // END READ DATA
150  // ==========================================================
151 
152 
153  _numPV = 0;
154  bool bPrimaryVertex = true;
155  if(doPVCheck_){
156  bPrimaryVertex = false;
157 
158  if (!vertex.isValid()) {
159  LogTrace(metname) << "[EfficiencyAnalyzer] Could not find vertex collection" << std::endl;
160  bPrimaryVertex = false;
161  }
162 
163  if ( vertex.isValid() ){
164  const reco::VertexCollection& vertexCollection = *(vertex.product());
165  int vertex_number = vertexCollection.size();
166 
167  reco::VertexCollection::const_iterator v = vertexCollection.begin();
168  for ( ; v != vertexCollection.end(); ++v) {
169  double vertex_chi2 = v->normalizedChi2();
170  double vertex_ndof = v->ndof();
171  bool fakeVtx = v->isFake();
172  double vertex_Z = v->z();
173 
174  if ( !fakeVtx
175  && vertex_number >= 1
176  && vertex_ndof > 4
177  && vertex_chi2 < 999
178  && fabs(vertex_Z)< 24. ) {
179  bPrimaryVertex = true;
180  ++_numPV;
181  }
182  }
183  }
184 
185  }
186 
187  // =================================================================================
188  // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
189  reco::Vertex::Point posVtx;
190  reco::Vertex::Error errVtx;
191  unsigned int theIndexOfThePrimaryVertex = 999.;
192  if (vertex.isValid()){
193  for (unsigned int ind=0; ind<vertex->size(); ++ind) {
194  if ( (*vertex)[ind].isValid() && !((*vertex)[ind].isFake()) ) {
195  theIndexOfThePrimaryVertex = ind;
196  break;
197  }
198  }
199  }
200 
201  if (theIndexOfThePrimaryVertex<100) {
202  posVtx = ((*vertex)[theIndexOfThePrimaryVertex]).position();
203  errVtx = ((*vertex)[theIndexOfThePrimaryVertex]).error();
204  }
205  else {
206  LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
207 
208  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
209  iEvent.getByToken(theBeamSpotLabel_,recoBeamSpotHandle);
210  reco::BeamSpot bs = *recoBeamSpotHandle;
211 
212  posVtx = bs.position();
213  errVtx(0,0) = bs.BeamWidthX();
214  errVtx(1,1) = bs.BeamWidthY();
215  errVtx(2,2) = bs.sigmaZ();
216  }
217 
218  const reco::Vertex thePrimaryVertex(posVtx,errVtx);
219  // ==========================================================
220 
221  if(!muons.isValid()) return;
222 
223  // Loop on muon collection
224  TLorentzVector Mu1, Mu2;
225 
226  bool isMB = false;
227  bool isME = false;
228 
229  for (reco::MuonCollection::const_iterator recoMu1 = muons->begin(); recoMu1!=muons->end(); ++recoMu1) {
230  LogTrace(metname)<<"[EfficiencyAnalyzer] loop over first muons" << endl;
231 
232  //--- Define combined isolation
233  reco::MuonIsolation Iso_muon = recoMu1->isolationR03();
234  float combIso = (Iso_muon.emEt + Iso_muon.hadEt + Iso_muon.sumPt);
235 
236  //--- Is Global Muon
237  if (!recoMu1->isGlobalMuon()) continue;
238 
239  // get the track combinig the information from both the Tracker and the Spectrometer
240  reco::TrackRef recoCombinedGlbTrack1 = recoMu1->combinedMuon();
241  float muPt1 = recoCombinedGlbTrack1->pt();
242  Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(), recoCombinedGlbTrack1->py(),recoCombinedGlbTrack1->pz(), recoCombinedGlbTrack1->p());
243 
244 
245  //--- Define if it is a tight muon
246  // Change the Tight muon definition by using the implemented method in: MuonSelectors.cc
247 
248  if (!muon::isTightMuon(*recoMu1, thePrimaryVertex)) continue;
249 
250  //-- is isolated muon
251  if (muPt1 <= 15) continue;
252  if (combIso/muPt1 > 0.1 ) continue;
253 
254  for (reco::MuonCollection::const_iterator recoMu2 = muons->begin(); recoMu2!=muons->end(); ++recoMu2){
255  LogTrace(metname)<<"[EfficiencyAnalyzer] loop over second muon" <<endl;
256  if (recoMu2 == recoMu1) continue;
257 
258  if (recoMu2->eta() < 1.479 ) isMB = true;
259  if (recoMu2->eta() >= 1.479 ) isME = true;
260 
261  //--> should we apply track quality cuts???
262  Mu2.SetPxPyPzE(recoMu2->px(), recoMu2->py(), recoMu2->pz(), recoMu2->p());
263 
264  float Minv = (Mu1+Mu2).M();
265  if (!recoMu2->isTrackerMuon()) continue;
266  if ( recoMu2->pt() < 5 ) continue;
267  if ( (recoMu1->charge())*(recoMu2->charge()) > 0 ) continue;
268  if ( Minv < 70 || Minv > 110 ) continue;
269 
270  h_allProbes_pt->Fill(recoMu2->pt());
271  h_allProbes_eta->Fill(recoMu2->eta());
272  h_allProbes_phi->Fill(recoMu2->phi());
273 
274  if (isMB) h_allProbes_EB_pt->Fill(recoMu2->pt());
275  if (isME) h_allProbes_EE_pt->Fill(recoMu2->pt());
276  if(recoMu2->pt() > 20 ) h_allProbes_hp_eta->Fill(recoMu2->eta());
277 
278  test_TightMu_Minv->Fill(Minv);
279 
280 
281 
282  // Probes passing the tight muon criteria
283 
284  if (!muon::isTightMuon(*recoMu2, thePrimaryVertex)) continue;
285 
286  h_passProbes_TightMu_pt->Fill(recoMu2->pt());
287  h_passProbes_TightMu_eta->Fill(recoMu2->eta());
288  h_passProbes_TightMu_phi->Fill(recoMu2->phi());
289 
290  if (isMB) h_passProbes_TightMu_EB_pt->Fill(recoMu2->pt());
291  if (isME) h_passProbes_TightMu_EE_pt->Fill(recoMu2->pt());
292  if( recoMu2->pt() > 20 ) h_passProbes_TightMu_hp_eta->Fill(recoMu2->eta());
293 
294  h_allProbes_TightMu_pt->Fill(recoMu2->pt());
295  if (isMB) h_allProbes_EB_TightMu_pt->Fill(recoMu2->pt());
296  if (isME) h_allProbes_EE_TightMu_pt->Fill(recoMu2->pt());
297 
298  //------- For PU monitoring -------//
299  if (bPrimaryVertex) h_allProbes_TightMu_nVtx->Fill(_numPV);
300  if (bPrimaryVertex && isMB) h_allProbes_EB_TightMu_nVtx->Fill(_numPV);
301  if (bPrimaryVertex && isME) h_allProbes_EE_TightMu_nVtx->Fill(_numPV);
302 
303  //-- Define det relative isolation
304  float tkIso = recoMu2->isolationR03().sumPt;
305  float emIso = recoMu2->isolationR03().emEt;
306  float hadIso = recoMu2->isolationR03().hadEt + recoMu2->isolationR03().hoEt;
307  float relDetIso = (tkIso + emIso + hadIso) / (recoMu2->pt());
308 
309  if (relDetIso < 0.05 ) {
310  h_passProbes_detIsoTightMu_pt->Fill(recoMu2->pt());
311  if (isMB) h_passProbes_EB_detIsoTightMu_pt->Fill(recoMu2->pt());
312  if (isME) h_passProbes_EE_detIsoTightMu_pt->Fill(recoMu2->pt());
313 
314  if (bPrimaryVertex) h_passProbes_detIsoTightMu_nVtx->Fill(_numPV);
315  if (bPrimaryVertex && isMB) h_passProbes_EB_detIsoTightMu_nVtx->Fill(_numPV);
316  if (bPrimaryVertex && isME) h_passProbes_EE_detIsoTightMu_nVtx->Fill(_numPV);
317  }
318 
319  //-- Define PF relative isolation
320  float chargedIso = recoMu2->pfIsolationR04().sumChargedHadronPt;
321  float neutralIso = recoMu2->pfIsolationR04().sumNeutralHadronEt;
322  float photonIso = recoMu2->pfIsolationR04().sumPhotonEt;
323  float relPFIso = (chargedIso + neutralIso + photonIso) / (recoMu2->pt());
324 
325  float pu = recoMu2->pfIsolationR04().sumPUPt;
326 
327  float neutralphotonPUCorrected = std::max(0.0 , (neutralIso + photonIso - 0.5*pu ) );
328 
329  float relPFIsoPUCorrected = (chargedIso + neutralphotonPUCorrected) / (recoMu2->pt());
330 
331 
332 
333  if (relPFIso < 0.12 ) {
334  h_passProbes_pfIsoTightMu_pt->Fill(recoMu2->pt());
335  if (isMB) h_passProbes_EB_pfIsoTightMu_pt->Fill(recoMu2->pt());
336  if (isME) h_passProbes_EE_pfIsoTightMu_pt->Fill(recoMu2->pt());
337 
338  if( bPrimaryVertex) h_passProbes_pfIsoTightMu_nVtx->Fill(_numPV);
339  if (bPrimaryVertex && isMB) h_passProbes_EB_pfIsoTightMu_nVtx->Fill(_numPV);
340  if (bPrimaryVertex && isME) h_passProbes_EE_pfIsoTightMu_nVtx->Fill(_numPV);
341  }
342 
343 
344 
345  // Apply deltaBeta PU corrections to the PF isolation eficiencies.
346 
347  if ( relPFIsoPUCorrected < 0.12 ) {
348 
349  h_passProbes_pfIsodBTightMu_pt->Fill(recoMu2->pt());
350  if (isMB) h_passProbes_EB_pfIsodBTightMu_pt->Fill(recoMu2->pt());
351  if (isME) h_passProbes_EE_pfIsodBTightMu_pt->Fill(recoMu2->pt());
352 
353  if( bPrimaryVertex) h_passProbes_pfIsodBTightMu_nVtx->Fill(_numPV);
354  if (bPrimaryVertex && isMB) h_passProbes_EB_pfIsodBTightMu_nVtx->Fill(_numPV);
355  if (bPrimaryVertex && isME) h_passProbes_EE_pfIsodBTightMu_nVtx->Fill(_numPV);
356  }
357  }
358  }
359 }
360 
361 
float hadEt
hcal sum-Et
Definition: MuonIsolation.h:9
dictionary parameters
Definition: Parameters.py:2
float sumPt
sum-pt of tracks
Definition: MuonIsolation.h:7
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
const std::string metname
math::Error< dimension >::type Error
covariance error matrix (3x3)
Definition: Vertex.h:43
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
tuple vertexCollection
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
int iEvent
Definition: GenABIO.cc:243
const T & max(const T &a, const T &b)
float emEt
ecal sum-Et
Definition: MuonIsolation.h:8
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:86
bool isValid() const
Definition: HandleBase.h:76
#define LogTrace(id)
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
void beginRun(const edm::Run &iRun, const edm::EventSetup &iSetup)
tuple tracks
Definition: testEve_cfg.py:39
double sigmaZ() const
sigma z
Definition: BeamSpot.h:80
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:88
T const * product() const
Definition: Handle.h:81
tuple muons
Definition: patZpeak.py:38
tuple cout
Definition: gather_cfg.py:121
const Point & position() const
position
Definition: BeamSpot.h:62
bool isTightMuon(const reco::Muon &, const reco::Vertex &)
EfficiencyAnalyzer(const edm::ParameterSet &pset)
Definition: Run.h:41