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