CMS 3D CMS Logo

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