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