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  ID_ = parameters.getParameter<string>("ID");
60 }
61 
63  delete theService;
64 }
65 
67  edm::Run const & /*iRun*/,
68  edm::EventSetup const & /* iSetup */){
69 
70  ibooker.cd();
71  ibooker.setCurrentFolder("Muons/EfficiencyAnalyzer/"+ID_);
72 
73  h_allProbes_pt = ibooker.book1D("allProbes_pt","All Probes Pt", ptBin_, ptMin_, ptMax_);
74  h_allProbes_EB_pt = ibooker.book1D("allProbes_EB_pt","Barrel: all Probes Pt", ptBin_, ptMin_, ptMax_);
75  h_allProbes_EE_pt = ibooker.book1D("allProbes_EE_pt","Endcap: all Probes Pt", ptBin_, ptMin_, ptMax_);
76  h_allProbes_eta = ibooker.book1D("allProbes_eta","All Probes Eta", etaBin_, etaMin_, etaMax_);
77  h_allProbes_hp_eta = ibooker.book1D("allProbes_hp_eta","High Pt all Probes Eta", etaBin_, etaMin_, etaMax_);
78  h_allProbes_phi = ibooker.book1D("allProbes_phi","All Probes Phi", phiBin_, phiMin_, phiMax_);
79 
80  h_allProbes_ID_pt = ibooker.book1D("allProbes_ID_pt","All ID Probes Pt", ptBin_, ptMin_, ptMax_);
81  h_allProbes_EB_ID_pt = ibooker.book1D("allProbes_EB_ID_pt","Barrel: all ID Probes Pt", ptBin_, ptMin_, ptMax_);
82  h_allProbes_EE_ID_pt = ibooker.book1D("allProbes_EE_ID_pt","Endcap: all ID Probes Pt", ptBin_, ptMin_, ptMax_);
83  h_allProbes_ID_nVtx = ibooker.book1D("allProbes_ID_nVtx","All Probes (ID) nVtx", vtxBin_, vtxMin_, vtxMax_);
84  h_allProbes_EB_ID_nVtx = ibooker.book1D("allProbes_EB_ID_nVtx","Barrel: All Probes (ID) nVtx", vtxBin_, vtxMin_, vtxMax_);
85  h_allProbes_EE_ID_nVtx = ibooker.book1D("allProbes_EE_ID_nVtx","Endcap: All Probes (ID) nVtx", vtxBin_, vtxMin_, vtxMax_);
86 
87  h_passProbes_ID_pt = ibooker.book1D("passProbes_ID_pt","ID Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
88  h_passProbes_ID_EB_pt = ibooker.book1D("passProbes_ID_EB_pt","Barrel: ID Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
89  h_passProbes_ID_EE_pt = ibooker.book1D("passProbes_ID_EE_pt","Endcap: ID Passing Probes Pt", ptBin_ , ptMin_ , ptMax_ );
90  h_passProbes_ID_eta = ibooker.book1D("passProbes_ID_eta","ID Passing Probes #eta", etaBin_, etaMin_, etaMax_);
91  h_passProbes_ID_hp_eta = ibooker.book1D("passProbes_ID_hp_eta","High Pt ID Passing Probes #eta", etaBin_, etaMin_, etaMax_);
92  h_passProbes_ID_phi = ibooker.book1D("passProbes_ID_phi","ID Passing Probes #phi", phiBin_, phiMin_, phiMax_);
93 
94  h_passProbes_detIsoID_pt = ibooker.book1D("passProbes_detIsoID_pt","detIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
95  h_passProbes_EB_detIsoID_pt = ibooker.book1D("passProbes_EB_detIsoID_pt","Barrel: detIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
96  h_passProbes_EE_detIsoID_pt = ibooker.book1D("passProbes_EE_detIsoID_pt","Endcap: detIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
97 
98  h_passProbes_pfIsoID_pt = ibooker.book1D("passProbes_pfIsoID_pt","pfIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
99  h_passProbes_EB_pfIsoID_pt = ibooker.book1D("passProbes_EB_pfIsoID_pt","Barrel: pfIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
100  h_passProbes_EE_pfIsoID_pt = ibooker.book1D("passProbes_EE_pfIsoID_pt","Endcap: pfIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
101 
102  h_passProbes_detIsoID_nVtx = ibooker.book1D("passProbes_detIsoID_nVtx", "detIsoID Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
103  h_passProbes_pfIsoID_nVtx = ibooker.book1D("passProbes_pfIsoID_nVtx", "pfIsoID Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
104  h_passProbes_EB_detIsoID_nVtx = ibooker.book1D("passProbes_EB_detIsoID_nVtx","Barrel: detIsoID Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
105  h_passProbes_EE_detIsoID_nVtx = ibooker.book1D("passProbes_EE_detIsoID_nVtx","Endcap: detIsoID Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
106  h_passProbes_EB_pfIsoID_nVtx = ibooker.book1D("passProbes_EB_pfIsoID_nVtx", "Barrel: pfIsoID Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
107  h_passProbes_EE_pfIsoID_nVtx = ibooker.book1D("passProbes_EE_pfIsoID_nVtx", "Endcap: pfIsoID Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
108 
109 
110  // Apply deltaBeta PU corrections to the PF isolation eficiencies.
111 
112  h_passProbes_pfIsodBID_pt = ibooker.book1D("passProbes_pfIsodBID_pt","pfIsoID Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
113  h_passProbes_EB_pfIsodBID_pt = ibooker.book1D("passProbes_EB_pfIsodBID_pt","Barrel: pfIsoID Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
114  h_passProbes_EE_pfIsodBID_pt = ibooker.book1D("passProbes_EE_pfIsodBID_pt","Endcap: pfIsoID Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
115  h_passProbes_pfIsodBID_nVtx = ibooker.book1D("passProbes_pfIsodBID_nVtx", "pfIsoID Passing Probes nVtx (R04) (deltaB PU correction)", vtxBin_, vtxMin_, vtxMax_);
116 h_passProbes_EB_pfIsodBID_nVtx = ibooker.book1D("passProbes_EB_pfIsodBID_nVtx", "Barrel: pfIsoID Passing Probes nVtx (R04) (deltaB PU correction)", vtxBin_, vtxMin_, vtxMax_);
117  h_passProbes_EE_pfIsodBID_nVtx = ibooker.book1D("passProbes_EE_pfIsodBID_nVtx", "Endcap: pfIsoID Passing Probes nVtx (R04) (deltaB PU correction)", vtxBin_, vtxMin_, vtxMax_);
118 
119 #ifdef DEBUG
120  cout << "[EfficiencyAnalyzer] Parameters initialization DONE" <<endl;
121 #endif
122 }
123 
125 
126  LogTrace(metname)<<"[EfficiencyAnalyzer] Analyze the mu in different eta regions";
127  theService->update(iSetup);
128  // ==========================================================
129  // BEGIN READ DATA:
130  // Muon information
132  iEvent.getByToken(theMuonCollectionLabel_, muons);
133 
134  // Tracks information
136  iEvent.getByToken(theTrackCollectionLabel_,tracks);
137 
138  //Vertex information
140  iEvent.getByToken(theVertexLabel_, vertex);
141  // END READ DATA
142  // ==========================================================
143 
144 
145  _numPV = 0;
146  bool bPrimaryVertex = true;
147  if(doPVCheck_){
148  bPrimaryVertex = false;
149 
150  if (!vertex.isValid()) {
151  LogTrace(metname) << "[EfficiencyAnalyzer] Could not find vertex collection" << std::endl;
152  bPrimaryVertex = false;
153  }
154 
155  if ( vertex.isValid() ){
156  const reco::VertexCollection& vertexCollection = *(vertex.product());
157  int vertex_number = vertexCollection.size();
158 
159  reco::VertexCollection::const_iterator v = vertexCollection.begin();
160  for ( ; v != vertexCollection.end(); ++v) {
161  double vertex_chi2 = v->normalizedChi2();
162  double vertex_ndof = v->ndof();
163  bool fakeVtx = v->isFake();
164  double vertex_Z = v->z();
165 
166  if ( !fakeVtx
167  && vertex_number >= 1
168  && vertex_ndof > 4
169  && vertex_chi2 < 999
170  && fabs(vertex_Z)< 24. ) {
171  bPrimaryVertex = true;
172  ++_numPV;
173  }
174  }
175  }
176 
177  }
178 
179  // =================================================================================
180  // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
181  reco::Vertex::Point posVtx;
182  reco::Vertex::Error errVtx;
183  unsigned int theIndexOfThePrimaryVertex = 999.;
184  if (vertex.isValid()){
185  for (unsigned int ind=0; ind<vertex->size(); ++ind) {
186  if ( (*vertex)[ind].isValid() && !((*vertex)[ind].isFake()) ) {
187  theIndexOfThePrimaryVertex = ind;
188  break;
189  }
190  }
191  }
192 
193  if (theIndexOfThePrimaryVertex<100) {
194  posVtx = ((*vertex)[theIndexOfThePrimaryVertex]).position();
195  errVtx = ((*vertex)[theIndexOfThePrimaryVertex]).error();
196  }
197  else {
198  LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
199 
200  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
201  iEvent.getByToken(theBeamSpotLabel_,recoBeamSpotHandle);
202  reco::BeamSpot bs = *recoBeamSpotHandle;
203 
204  posVtx = bs.position();
205  errVtx(0,0) = bs.BeamWidthX();
206  errVtx(1,1) = bs.BeamWidthY();
207  errVtx(2,2) = bs.sigmaZ();
208  }
209 
210  const reco::Vertex thePrimaryVertex(posVtx,errVtx);
211  // ==========================================================
212 
213  if(!muons.isValid()) return;
214 
215  // Loop on muon collection
216  TLorentzVector Mu1, Mu2;
217 
218  bool isMB = false;
219  bool isME = false;
220 
221  for (reco::MuonCollection::const_iterator recoMu1 = muons->begin(); recoMu1!=muons->end(); ++recoMu1) {
222  LogTrace(metname)<<"[EfficiencyAnalyzer] loop over first muons" << endl;
223 
224  //--- Define combined isolation
225  reco::MuonIsolation Iso_muon = recoMu1->isolationR03();
226  float combIso = (Iso_muon.emEt + Iso_muon.hadEt + Iso_muon.sumPt);
227 
228  //--- Is Global Muon
229  if (!recoMu1->isGlobalMuon()) continue;
230 
231  // get the track combinig the information from both the Tracker and the Spectrometer
232  reco::TrackRef recoCombinedGlbTrack1 = recoMu1->combinedMuon();
233  float muPt1 = recoCombinedGlbTrack1->pt();
234  Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(), recoCombinedGlbTrack1->py(),recoCombinedGlbTrack1->pz(), recoCombinedGlbTrack1->p());
235 
236 
237  //--- Define if it is a tight muon
238  // Change the Tight muon definition by using the implemented method in: MuonSelectors.cc
239  if (ID_ == "Loose" && !muon::isLooseMuon(*recoMu1)) continue;
240  if (ID_ == "Medium" && !muon::isMediumMuon(*recoMu1)) continue;
241  if (ID_ == "Tight" && !muon::isTightMuon(*recoMu1, thePrimaryVertex)) continue;
242 
243  //-- is isolated muon
244  if (muPt1 <= 15) continue;
245  if (combIso/muPt1 > 0.1 ) continue;
246 
247  for (reco::MuonCollection::const_iterator recoMu2 = muons->begin(); recoMu2!=muons->end(); ++recoMu2){
248  LogTrace(metname)<<"[EfficiencyAnalyzer] loop over second muon" <<endl;
249  if (recoMu2 == recoMu1) continue;
250 
251  if (recoMu2->eta() < 1.479 ) isMB = true;
252  if (recoMu2->eta() >= 1.479 ) isME = true;
253 
254  //--> should we apply track quality cuts???
255  Mu2.SetPxPyPzE(recoMu2->px(), recoMu2->py(), recoMu2->pz(), recoMu2->p());
256 
257  float Minv = (Mu1+Mu2).M();
258  if (!recoMu2->isTrackerMuon()) continue;
259  if ( recoMu2->pt() < 5 ) continue;
260  if ( (recoMu1->charge())*(recoMu2->charge()) > 0 ) continue;
261  if ( Minv < 70 || Minv > 110 ) continue;
262 
263  h_allProbes_pt->Fill(recoMu2->pt());
264  h_allProbes_eta->Fill(recoMu2->eta());
265  h_allProbes_phi->Fill(recoMu2->phi());
266 
267  if (isMB) h_allProbes_EB_pt->Fill(recoMu2->pt());
268  if (isME) h_allProbes_EE_pt->Fill(recoMu2->pt());
269  if(recoMu2->pt() > 20 ) h_allProbes_hp_eta->Fill(recoMu2->eta());
270 
271 
272  // Probes passing the tight muon criteria
273  if (ID_ == "Loose" && !muon::isLooseMuon(*recoMu2)) continue;
274  if (ID_ == "Medium" && !muon::isMediumMuon(*recoMu2)) continue;
275  if (ID_ == "Tight" && !muon::isTightMuon(*recoMu2, thePrimaryVertex)) continue;
276 
277 
278  h_passProbes_ID_pt->Fill(recoMu2->pt());
279  h_passProbes_ID_eta->Fill(recoMu2->eta());
280  h_passProbes_ID_phi->Fill(recoMu2->phi());
281 
282  if (isMB) h_passProbes_ID_EB_pt->Fill(recoMu2->pt());
283  if (isME) h_passProbes_ID_EE_pt->Fill(recoMu2->pt());
284  if( recoMu2->pt() > 20 ) h_passProbes_ID_hp_eta->Fill(recoMu2->eta());
285 
286  h_allProbes_ID_pt->Fill(recoMu2->pt());
287  if (isMB) h_allProbes_EB_ID_pt->Fill(recoMu2->pt());
288  if (isME) h_allProbes_EE_ID_pt->Fill(recoMu2->pt());
289 
290  //------- For PU monitoring -------//
291  if (bPrimaryVertex) h_allProbes_ID_nVtx->Fill(_numPV);
292  if (bPrimaryVertex && isMB) h_allProbes_EB_ID_nVtx->Fill(_numPV);
293  if (bPrimaryVertex && isME) h_allProbes_EE_ID_nVtx->Fill(_numPV);
294 
295  //-- Define det relative isolation
296  float tkIso = recoMu2->isolationR03().sumPt;
297  float emIso = recoMu2->isolationR03().emEt;
298  float hadIso = recoMu2->isolationR03().hadEt + recoMu2->isolationR03().hoEt;
299  float relDetIso = (tkIso + emIso + hadIso) / (recoMu2->pt());
300 
301  if (relDetIso < 0.05 ) {
302  h_passProbes_detIsoID_pt->Fill(recoMu2->pt());
303  if (isMB) h_passProbes_EB_detIsoID_pt->Fill(recoMu2->pt());
304  if (isME) h_passProbes_EE_detIsoID_pt->Fill(recoMu2->pt());
305 
306  if (bPrimaryVertex) h_passProbes_detIsoID_nVtx->Fill(_numPV);
307  if (bPrimaryVertex && isMB) h_passProbes_EB_detIsoID_nVtx->Fill(_numPV);
308  if (bPrimaryVertex && isME) h_passProbes_EE_detIsoID_nVtx->Fill(_numPV);
309  }
310 
311  //-- Define PF relative isolation
312  float chargedIso = recoMu2->pfIsolationR04().sumChargedHadronPt;
313  float neutralIso = recoMu2->pfIsolationR04().sumNeutralHadronEt;
314  float photonIso = recoMu2->pfIsolationR04().sumPhotonEt;
315  float relPFIso = (chargedIso + neutralIso + photonIso) / (recoMu2->pt());
316 
317  float pu = recoMu2->pfIsolationR04().sumPUPt;
318  float neutralphotonPUCorrected = std::max(0.0 , (neutralIso + photonIso - 0.5*pu ) );
319  float relPFIsoPUCorrected = (chargedIso + neutralphotonPUCorrected) / (recoMu2->pt());
320 
321 
322 
323  if (relPFIso < 0.12 ) {
324  h_passProbes_pfIsoID_pt->Fill(recoMu2->pt());
325  if (isMB) h_passProbes_EB_pfIsoID_pt->Fill(recoMu2->pt());
326  if (isME) h_passProbes_EE_pfIsoID_pt->Fill(recoMu2->pt());
327 
328  if( bPrimaryVertex) h_passProbes_pfIsoID_nVtx->Fill(_numPV);
329  if (bPrimaryVertex && isMB) h_passProbes_EB_pfIsoID_nVtx->Fill(_numPV);
330  if (bPrimaryVertex && isME) h_passProbes_EE_pfIsoID_nVtx->Fill(_numPV);
331  }
332 
333  // Apply deltaBeta PU corrections to the PF isolation eficiencies.
334  if ( relPFIsoPUCorrected < 0.12 ) {
335  h_passProbes_pfIsodBID_pt->Fill(recoMu2->pt());
336  if (isMB) h_passProbes_EB_pfIsodBID_pt->Fill(recoMu2->pt());
337  if (isME) h_passProbes_EE_pfIsodBID_pt->Fill(recoMu2->pt());
338 
339  if( bPrimaryVertex) h_passProbes_pfIsodBID_nVtx->Fill(_numPV);
340  if (bPrimaryVertex && isMB) h_passProbes_EB_pfIsodBID_nVtx->Fill(_numPV);
341  if (bPrimaryVertex && isME) h_passProbes_EE_pfIsodBID_nVtx->Fill(_numPV);
342  }
343  }
344  }
345 }
346 
347 
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:449
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
bool isMediumMuon(const reco::Muon &)
tuple vertexCollection
bool isLooseMuon(const reco::Muon &)
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
int iEvent
Definition: GenABIO.cc:230
float emEt
ecal sum-Et
Definition: MuonIsolation.h:8
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
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:75
#define LogTrace(id)
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
T const * product() const
Definition: Handle.h:81
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
tuple muons
Definition: patZpeak.py:38
static int position[264][3]
Definition: ReadPGInfo.cc:509
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