CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ZMuMuEfficiency.cc
Go to the documentation of this file.
1 /* \class ZMuMuEfficiency
2  *
3  * author: Pasquale Noli
4  * revised by Salvatore di Guida
5  * revised for CSA08 by Davide Piccolo
6  *
7  * Efficiency of reconstruction tracker and muon Chamber
8  *
9  */
10 
22 #include "TH1.h"
23 #include <vector>
24 using namespace edm;
25 using namespace std;
26 using namespace reco;
27 
29 public:
31 private:
32  virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
33  bool check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
34  float getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
35  float getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
36  float getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
37  Particle::LorentzVector getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
38  virtual void endJob() override;
39 
43  edm::InputTag muons_, muonMatchMap_, muonIso_;
44  edm::InputTag tracks_, trackIso_;
47 
48  double zMassMin_, zMassMax_, ptmin_, etamax_, isomax_;
49  unsigned int nbinsPt_, nbinsEta_;
50  reco::CandidateBaseRef globalMuonCandRef_, trackMuonCandRef_, standAloneMuonCandRef_;
52 
53  //histograms for measuring tracker efficiency
54  TH1D *h_etaStandAlone_, *h_etaMuonOverlappedToStandAlone_;
55  TH1D *h_ptStandAlone_, *h_ptMuonOverlappedToStandAlone_;
56 
57  //histograms for measuring standalone efficiency
58  TH1D *h_etaTrack_, *h_etaMuonOverlappedToTrack_;
59  TH1D *h_ptTrack_, *h_ptMuonOverlappedToTrack_;
60 
61  //histograms for MC acceptance
62  TH1D *h_nZMCfound_;
63  TH1D *h_ZetaGen_, *h_ZptGen_, *h_ZmassGen_;
64  TH1D *h_muetaGen_, *h_muptGen_, *h_muIsoGen_;
65  TH1D *h_dimuonPtGen_, *h_dimuonMassGen_, *h_dimuonEtaGen_;
66  TH1D *h_ZetaGenPassed_, *h_ZptGenPassed_, *h_ZmassGenPassed_;
67  TH1D *h_muetaGenPassed_, *h_muptGenPassed_, *h_muIsoGenPassed_;
68  TH1D *h_dimuonPtGenPassed_, *h_dimuonMassGenPassed_, *h_dimuonEtaGenPassed_;
69  //histograms for invarian mass resolution
70  TH1D *h_DELTA_ZMuMuMassReco_dimuonMassGen_, *h_DELTA_ZMuStaMassReco_dimuonMassGen_, *h_DELTA_ZMuTrackMassReco_dimuonMassGen_;
71 
72  int numberOfEventsWithZMuMufound, numberOfEventsWithZMuStafound;
73  int numberOfMatchedZMuSta_,numberOfMatchedSelectedZMuSta_;
74  int numberOfMatchedZMuMu_, numberOfMatchedSelectedZMuMu_;
75  int numberOfOverlappedStandAlone_, numberOfOverlappedTracks_, numberOfMatchedZMuTrack_notOverlapped;
76  int numberOfMatchedZMuTrack_exclusive ,numberOfMatchedSelectedZMuTrack_exclusive;
77  int numberOfMatchedZMuTrack_matchedZMuMu, numberOfMatchedZMuTrack_matchedSelectedZMuMu ;
78  int totalNumberOfevents, totalNumberOfZfound, totalNumberOfZPassed;
79  int noMCmatching, ZMuTrack_exclusive_1match, ZMuTrack_exclusive_morematch;
80  int ZMuTrackselected_exclusive_1match, ZMuTrackselected_exclusive_morematch;
81  int ZMuTrack_ZMuMu_1match, ZMuTrack_ZMuMu_2match, ZMuTrack_ZMuMu_morematch;
82 
83  int n_zMuMufound_genZsele, n_zMuStafound_genZsele, n_zMuTrkfound_genZsele;
84 };
85 
94 #include <iostream>
95 #include <iterator>
96 #include <cmath>
97 using namespace std;
98 using namespace reco;
99 using namespace edm;
100 
101 
103 
105  zMuMu_(pset.getParameter<InputTag>("zMuMu")),
106  zMuMuMatchMap_(pset.getParameter<InputTag>("zMuMuMatchMap")),
107  zMuTrack_(pset.getParameter<InputTag>("zMuTrack")),
108  zMuTrackMatchMap_(pset.getParameter<InputTag>("zMuTrackMatchMap")),
109  zMuStandAlone_(pset.getParameter<InputTag>("zMuStandAlone")),
110  zMuStandAloneMatchMap_(pset.getParameter<InputTag>("zMuStandAloneMatchMap")),
111  muons_(pset.getParameter<InputTag>("muons")),
112  muonMatchMap_(pset.getParameter<InputTag>("muonMatchMap")),
113  muonIso_(pset.getParameter<InputTag>("muonIso")),
114  tracks_(pset.getParameter<InputTag>("tracks")),
115  trackIso_(pset.getParameter<InputTag>("trackIso")),
116  standAlone_(pset.getParameter<InputTag>("standAlone")),
117  standAloneIso_(pset.getParameter<InputTag>("standAloneIso")),
118  genParticles_(pset.getParameter<InputTag>( "genParticles" ) ),
119 
120  zMassMin_(pset.getUntrackedParameter<double>("zMassMin")),
121  zMassMax_(pset.getUntrackedParameter<double>("zMassMax")),
122  ptmin_(pset.getUntrackedParameter<double>("ptmin")),
123  etamax_(pset.getUntrackedParameter<double>("etamax")),
124  isomax_(pset.getUntrackedParameter<double>("isomax")),
125  nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
126  nbinsEta_(pset.getUntrackedParameter<unsigned int>("nbinsEta")) {
128  TFileDirectory trackEffDir = fs->mkdir("TrackEfficiency");
129 
130  // tracker efficiency distributions
131  h_etaStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonEta",
132  "StandAlone #eta for Z -> #mu + standalone",
134  h_etaMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAloneEta",
135  "Global muon overlapped to standAlone #eta for Z -> #mu + sa",
137  h_ptStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonPt",
138  "StandAlone p_{t} for Z -> #mu + standalone",
139  nbinsPt_, ptmin_, 100);
140  h_ptMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAlonePt",
141  "Global muon overlapped to standAlone p_{t} for Z -> #mu + sa",
142  nbinsPt_, ptmin_, 100);
143 
144 
145  // StandAlone efficiency distributions
146  TFileDirectory standaloneEffDir = fs->mkdir("StandaloneEfficiency");
147  h_etaTrack_ = standaloneEffDir.make<TH1D>("TrackMuonEta",
148  "Track #eta for Z -> #mu + track",
150  h_etaMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackEta",
151  "Global muon overlapped to track #eta for Z -> #mu + tk",
153  h_ptTrack_ = standaloneEffDir.make<TH1D>("TrackMuonPt",
154  "Track p_{t} for Z -> #mu + track",
155  nbinsPt_, ptmin_, 100);
156  h_ptMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackPt",
157  "Global muon overlapped to track p_{t} for Z -> #mu + tk",
158  nbinsPt_, ptmin_, 100);
159 
160 
161  // inv. mass resolution studies
162  TFileDirectory invMassResolutionDir = fs->mkdir("invriantMassResolution");
163  h_DELTA_ZMuMuMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuMu_invMassResolution","zMuMu invariant Mass Resolution",50,-25,25);
164  h_DELTA_ZMuStaMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuSta_invMassResolution","zMuSta invariant Mass Resolution",50,-25,25);
165  h_DELTA_ZMuTrackMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuTrack_invMassResolution","zMuTrack invariant Mass Resolution",50,-25,25);
166 
167 
168  // generator level histograms
169  TFileDirectory genParticleDir = fs->mkdir("genParticle");
170  h_nZMCfound_ = genParticleDir.make<TH1D>("NumberOfgeneratedZeta","n. of generated Z per event",4,-.5,3.5);
171  h_ZetaGen_ = genParticleDir.make<TH1D>("generatedZeta","#eta of generated Z",100,-5.,5.);
172  h_ZptGen_ = genParticleDir.make<TH1D>("generatedZpt","pt of generated Z",100,0.,200.);
173  h_ZmassGen_ = genParticleDir.make<TH1D>("generatedZmass","mass of generated Z",100,0.,200.);
174  h_muetaGen_ = genParticleDir.make<TH1D>("generatedMuonEta","#eta of generated muons from Z decay",100,-5.,5.);
175  h_muptGen_ = genParticleDir.make<TH1D>("generatedMuonpt","pt of generated muons from Z decay",100,0.,200.);
176  h_dimuonEtaGen_ = genParticleDir.make<TH1D>("generatedDimuonEta","#eta of generated dimuon",100,-5.,5.);
177  h_dimuonPtGen_ = genParticleDir.make<TH1D>("generatedDimuonPt","pt of generated dimuon",100,0.,200.);
178  h_dimuonMassGen_ = genParticleDir.make<TH1D>("generatedDimuonMass","mass of generated dimuon",100,0.,200.);
179  h_ZetaGenPassed_ = genParticleDir.make<TH1D>("generatedZeta_passed","#eta of generated Z after cuts",100,-5.,5.);
180  h_ZptGenPassed_ = genParticleDir.make<TH1D>("generatedZpt_passed","pt of generated Z after cuts",100,0.,200.);
181  h_ZmassGenPassed_ = genParticleDir.make<TH1D>("generatedZmass_passed","mass of generated Z after cuts",100,0.,200.);
182  h_muetaGenPassed_ = genParticleDir.make<TH1D>("generatedMuonEta_passed","#eta of generated muons from Z decay after cuts",100,-5.,5.);
183  h_muptGenPassed_ = genParticleDir.make<TH1D>("generatedMuonpt_passed","pt of generated muons from Z decay after cuts",100,0.,200.);
184  h_dimuonEtaGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonEta_passed","#eta of generated dimuon after cuts",100,-5.,5.);
185  h_dimuonPtGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonPt_passed","pt of generated dimuon after cuts",100,0.,200.);
186  h_dimuonMassGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonMass_passed","mass of generated dimuon after cuts",100,0.,200.);
187  // to insert isolation histograms ..............
188 
202  noMCmatching = 0;
210 
214 
215  // generator counters
219 
220 }
221 
224  Handle<GenParticleMatch> zMuMuMatchMap; //Map of Z made by Mu global + Mu global
225  Handle<CandidateView> zMuTrack;
226  Handle<GenParticleMatch> zMuTrackMatchMap; //Map of Z made by Mu + Track
227  Handle<CandidateView> zMuStandAlone;
228  Handle<GenParticleMatch> zMuStandAloneMatchMap; //Map of Z made by Mu + StandAlone
229  Handle<CandidateView> muons; //Collection of Muons
230  Handle<GenParticleMatch> muonMatchMap;
232  Handle<CandidateView> tracks; //Collection of Tracks
233  Handle<IsolationCollection> trackIso;
234  Handle<CandidateView> standAlone; //Collection of StandAlone
235  Handle<IsolationCollection> standAloneIso;
236  Handle<GenParticleCollection> genParticles; // Collection of Generatd Particles
237 
238  event.getByLabel(zMuMu_, zMuMu);
239  event.getByLabel(zMuTrack_, zMuTrack);
240  event.getByLabel(zMuStandAlone_, zMuStandAlone);
241  event.getByLabel(muons_, muons);
242  event.getByLabel(tracks_, tracks);
243  event.getByLabel(standAlone_, standAlone);
244  event.getByLabel(genParticles_, genParticles);
245 
246  cout << "********* zMuMu size : " << zMuMu->size() << endl;
247  cout << "********* zMuStandAlone size : " << zMuStandAlone->size() << endl;
248  cout << "********* zMuTrack size : " << zMuTrack->size() << endl;
249  cout << "********* muons size : " << muons->size()<< endl;
250  cout << "********* standAlone size : " << standAlone->size()<< endl;
251  cout << "********* tracks size : " << tracks->size()<< endl;
252  cout << "********* generated size : " << genParticles->size()<< endl;
253 
254 
255  // generator level distributions
256 
257  int nZMCfound = 0;
259  int ngen = genParticles->size();
260  bool ZMuMuMatchedfound = false;
261  bool ZMuMuMatchedSelectedfound = false;
262  bool ZMuStaMatchedfound = false;
263  //bool ZMuStaMatchedSelectedfound = false;
264  int ZMuTrackMatchedfound = 0;
265  int ZMuTrackMatchedSelected_exclusivefound = 0;
266 
267  double dimuonMassGen = 0;
268 
269  for (int i=0; i<ngen; i++) {
270  const Candidate &genCand = (*genParticles)[i];
271 
272  // if((genCand.pdgId() == 23) && (genCand.status() == 2)) //this is an intermediate Z0
273  // cout << ">>> intermediate Z0 found, with " << genCand.numberOfDaughters() << " daughters" << endl;
274  if((genCand.pdgId() == 23)&&(genCand.status() == 3)) { //this is a Z0
275  if(genCand.numberOfDaughters() == 3) { // possible Z0 decays in mu+ mu-, the 3rd daughter is the same Z0
276  const Candidate * dauGen0 = genCand.daughter(0);
277  const Candidate * dauGen1 = genCand.daughter(1);
278  const Candidate * dauGen2 = genCand.daughter(2);
279  if (check_ifZmumu(dauGen0, dauGen1, dauGen2)) { // Z0 in mu+ mu-
281  nZMCfound++;
282  bool checkpt = false;
283  bool checketa = false;
284  bool checkmass = false;
285  float mupluspt, muminuspt, mupluseta, muminuseta;
286  mupluspt = getParticlePt(-13,dauGen0,dauGen1,dauGen2);
287  muminuspt = getParticlePt(13,dauGen0,dauGen1,dauGen2);
288  mupluseta = getParticleEta(-13,dauGen0,dauGen1,dauGen2);
289  muminuseta = getParticleEta(13,dauGen0,dauGen1,dauGen2);
290  //float muplusphi = getParticlePhi(-13,dauGen0,dauGen1,dauGen2);
291  //float muminusphi = getParticlePhi(13,dauGen0,dauGen1,dauGen2);
292  Particle::LorentzVector pZ(0, 0, 0, 0);
293  Particle::LorentzVector muplusp4 = getParticleP4(-13,dauGen0,dauGen1,dauGen2);
294  Particle::LorentzVector muminusp4 = getParticleP4(13,dauGen0,dauGen1,dauGen2);
295  pZ = muplusp4 + muminusp4;
296  double dimuon_pt = sqrt(pZ.x()*pZ.x()+pZ.y()*pZ.y());
297  double tan_theta_half = tan(atan(dimuon_pt/pZ.z())/2.);
298  double dimuon_eta = 0.;
299  if (tan_theta_half>0) dimuon_eta = -log(tan(tan_theta_half));
300  if (tan_theta_half<=0) dimuon_eta = log(tan(-tan_theta_half));
301 
302  dimuonMassGen = pZ.mass(); // dimuon invariant Mass at Generator Level
303 
304  h_ZmassGen_->Fill(genCand.mass());
305  h_ZetaGen_->Fill(genCand.eta());
306  h_ZptGen_->Fill(genCand.pt());
307  h_dimuonMassGen_->Fill(pZ.mass());
308  h_dimuonEtaGen_->Fill(dimuon_eta);
309  h_dimuonPtGen_->Fill(dimuon_pt);
310  h_muetaGen_->Fill(mupluseta);
311  h_muetaGen_->Fill(muminuseta);
312  h_muptGen_->Fill(mupluspt);
313  h_muptGen_->Fill(muminuspt);
314  // dimuon 4-momentum
315  // h_mDimuonMC->Fill(pZ.mass());
316  // h_ZminusDimuonMassMC->Fill(genCand.mass()-pZ.mass());
317  // h_DeltaPhiMC->Fill(deltaPhi(muplusphi,muminusphi));
318  // if (dauGen2==23) float z_eta = dauGen2->eta();
319  // if (dauGen2==23) float Zpt = dauGen2->pt();
320  // h_DeltaPhivsZPtMC->Fill(DeltaPhi(muplusphi,muminusphi),ZPt);
321 
322  if (mupluspt > ptmin_ && muminuspt > ptmin_) checkpt = true;
323  if (mupluseta < etamax_ && muminuseta < etamax_) checketa = true;
324  if (genCand.mass()>zMassMin_ && genCand.mass()<zMassMax_) checkmass = true;
325  if (checkpt && checketa && checkmass) {
327  h_ZmassGenPassed_->Fill(genCand.mass());
328  h_ZetaGenPassed_->Fill(genCand.eta());
329  h_ZptGenPassed_->Fill(genCand.pt());
330  h_dimuonMassGenPassed_->Fill(pZ.mass());
331  h_dimuonEtaGenPassed_->Fill(dimuon_eta);
332  h_dimuonPtGenPassed_->Fill(dimuon_pt);
333  h_muetaGenPassed_->Fill(mupluseta);
334  h_muetaGenPassed_->Fill(muminuseta);
335  h_muptGenPassed_->Fill(mupluspt);
336  h_muptGenPassed_->Fill(muminuspt);
337 
338  if (zMuMu->size() > 0 ) {
340  }
341  else if (zMuStandAlone->size() > 0 ) {
343  }
344  else {
346  }
347 
348  }
349  }
350 
351  }
352  }
353  }
354  h_nZMCfound_->Fill(nZMCfound); // number of Z found in the event at generator level
355 
356  //TRACK efficiency (conto numero di eventi Zmumu global e ZmuSta (ricorda che sono due campioni esclusivi)
357 
358  if (zMuMu->size() > 0 ) {
360  event.getByLabel(zMuMuMatchMap_, zMuMuMatchMap);
361  event.getByLabel(muonIso_, muonIso);
362  event.getByLabel(standAloneIso_, standAloneIso);
363  event.getByLabel(muonMatchMap_, muonMatchMap);
364  for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates
365  const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate
366  CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i);
367  bool isMatched = false;
368  GenParticleRef zMuMuMatch = (*zMuMuMatchMap)[zMuMuCandRef];
369 
370  if(zMuMuMatch.isNonnull()) { // ZMuMu matched
371  isMatched = true;
373  }
374  CandidateBaseRef dau0 = zMuMuCand.daughter(0)->masterClone();
375  CandidateBaseRef dau1 = zMuMuCand.daughter(1)->masterClone();
376  if (isMatched) ZMuMuMatchedfound = true;
377 
378  // Cuts
379  if((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) &&
380  (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta()) < etamax_) &&
381  (zMuMuCand.mass() > zMassMin_) && (zMuMuCand.mass() < zMassMax_) &&
382  (isMatched)) {
383  //The Z daughters are already matched!
384  const double globalMuonIsolation0 = (*muonIso)[dau0];
385  const double globalMuonIsolation1 = (*muonIso)[dau1];
386  if((globalMuonIsolation0 < isomax_) && (globalMuonIsolation1 < isomax_)) { // ZMuMu matched and selected by cuts
387  ZMuMuMatchedSelectedfound = true;
389  h_etaStandAlone_->Fill(dau0->eta()); // StandAlone found dau0, eta
390  h_etaStandAlone_->Fill(dau1->eta()); // StandAlone found dau1, eta
391  h_etaMuonOverlappedToStandAlone_->Fill(dau0->eta()); // is global muon so dau0 is also found as a track, eta
392  h_etaMuonOverlappedToStandAlone_->Fill(dau1->eta()); // is global muon so dau1 is also found as a track, eta
393  h_ptStandAlone_->Fill(dau0->pt()); // StandAlone found dau0, pt
394  h_ptStandAlone_->Fill(dau1->pt()); // StandAlone found dau1, pt
395  h_ptMuonOverlappedToStandAlone_->Fill(dau0->pt()); // is global muon so dau0 is also found as a track, pt
396  h_ptMuonOverlappedToStandAlone_->Fill(dau1->pt()); // is global muon so dau1 is also found as a track, pt
397 
398  h_etaTrack_->Fill(dau0->eta()); // Track found dau0, eta
399  h_etaTrack_->Fill(dau1->eta()); // Track found dau1, eta
400  h_etaMuonOverlappedToTrack_->Fill(dau0->eta()); // is global muon so dau0 is also found as a StandAlone, eta
401  h_etaMuonOverlappedToTrack_->Fill(dau1->eta()); // is global muon so dau1 is also found as a StandAlone, eta
402  h_ptTrack_->Fill(dau0->pt()); // Track found dau0, pt
403  h_ptTrack_->Fill(dau1->pt()); // Track found dau1, pt
404  h_ptMuonOverlappedToTrack_->Fill(dau0->pt()); // is global muon so dau0 is also found as a StandAlone, pt
405  h_ptMuonOverlappedToTrack_->Fill(dau1->pt()); // is global muon so dau1 is also found as a StandAlone, pt
406 
407  h_DELTA_ZMuMuMassReco_dimuonMassGen_->Fill(zMuMuCand.mass()-dimuonMassGen);
408  // check that the two muons are matched . .per ora è solo un mio controllo
409  for(unsigned int j = 0; j < muons->size() ; ++j) {
410  CandidateBaseRef muCandRef = muons->refAt(j);
411  GenParticleRef muonMatch = (*muonMatchMap)[muCandRef];
412  // if (muonMatch.isNonnull()) cout << "mu match n. " << j << endl;
413  }
414  }
415  }
416  }
417  }
418 
419  if (zMuStandAlone->size() > 0) {
421  event.getByLabel(zMuStandAloneMatchMap_, zMuStandAloneMatchMap);
422  event.getByLabel(muonIso_, muonIso);
423  event.getByLabel(standAloneIso_, standAloneIso);
424  event.getByLabel(muonMatchMap_, muonMatchMap);
425  for(unsigned int i = 0; i < zMuStandAlone->size(); ++i) { //loop on candidates
426  const Candidate & zMuStaCand = (*zMuStandAlone)[i]; //the candidate
427  CandidateBaseRef zMuStaCandRef = zMuStandAlone->refAt(i);
428  bool isMatched = false;
429  GenParticleRef zMuStaMatch = (*zMuStandAloneMatchMap)[zMuStaCandRef];
430  if(zMuStaMatch.isNonnull()) { // ZMuSta Macthed
431  isMatched = true;
432  ZMuStaMatchedfound = true;
434  }
435  CandidateBaseRef dau0 = zMuStaCand.daughter(0)->masterClone();
436  CandidateBaseRef dau1 = zMuStaCand.daughter(1)->masterClone();
437 
438  // Cuts
439  if((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) &&
440  (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta()) < etamax_) &&
441  (zMuStaCand.mass() > zMassMin_) && (zMuStaCand.mass() < zMassMax_) &&
442  (isMatched)) {
444  if(dau0->isGlobalMuon()) {
445  standAloneMuonCandRef_ = dau1;
446  globalMuonCandRef_ = dau0;
447  }
448  if(dau1->isGlobalMuon()) {
449  standAloneMuonCandRef_ = dau0;
450  globalMuonCandRef_ = dau1;
451  }
452 
453  const double globalMuonIsolation = (*muonIso)[globalMuonCandRef_];
454  const double standAloneMuonIsolation = (*standAloneIso)[standAloneMuonCandRef_];
455 
456  if((globalMuonIsolation < isomax_) && (standAloneMuonIsolation < isomax_)) { // ZMuSta matched and selected
457  //ZMuStaMatchedSelectedfound = true;
459  h_etaStandAlone_->Fill(standAloneMuonCandRef_->eta()); //Denominator eta for measuring track efficiency
460  h_ptStandAlone_->Fill(standAloneMuonCandRef_->pt()); //Denominator pt for measuring track eff
461  h_DELTA_ZMuStaMassReco_dimuonMassGen_->Fill(zMuStaCand.mass()-dimuonMassGen); // differnce between ZMuSta reco and dimuon mass gen
462 
463  }
464  }
465  }
466  } //end loop on Candidate
467 
468  //STANDALONE efficiency
469 
470  if (zMuTrack->size() > 0) {
471  event.getByLabel(zMuTrackMatchMap_, zMuTrackMatchMap);
472  event.getByLabel(muonIso_, muonIso);
473  event.getByLabel(trackIso_, trackIso);
474  event.getByLabel(muonMatchMap_, muonMatchMap);
475  for(unsigned int i = 0; i < zMuTrack->size(); ++i) { //loop on candidates
476  const Candidate & zMuTrkCand = (*zMuTrack)[i]; //the candidate
477  CandidateBaseRef zMuTrkCandRef = zMuTrack->refAt(i);
478  bool isMatched = false;
479  GenParticleRef zMuTrkMatch = (*zMuTrackMatchMap)[zMuTrkCandRef];
480  if(zMuTrkMatch.isNonnull()) {
481  isMatched = true;
482  }
483  CandidateBaseRef dau0 = zMuTrkCand.daughter(0)->masterClone();
484  CandidateBaseRef dau1 = zMuTrkCand.daughter(1)->masterClone();
485 
486  if (isMatched) {
487  ZMuTrackMatchedfound++;
488  if (ZMuMuMatchedfound) numberOfMatchedZMuTrack_matchedZMuMu++;
489  if (ZMuMuMatchedSelectedfound) numberOfMatchedZMuTrack_matchedSelectedZMuMu++;
490  if (!ZMuMuMatchedfound) numberOfMatchedZMuTrack_exclusive++;
491  }
492  // Cuts
493  if ((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) &&
494  (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta())< etamax_) &&
495  (zMuTrkCand.mass() > zMassMin_) && (zMuTrkCand.mass() < zMassMax_) &&
496  (isMatched) && !ZMuMuMatchedfound && !ZMuStaMatchedfound ) {
497 
498  // dau0 is always the global muon, dau1 is the track for ZMuTrack collection
499  const double globalMuonIsolation = (*muonIso)[dau0];
500  const double trackMuonIsolation = (*trackIso)[dau1];
501  if((globalMuonIsolation < isomax_) && (trackMuonIsolation < isomax_)) { // ZMuTRack matched - selected without ZMuMu found (exclusive)
503  ZMuTrackMatchedSelected_exclusivefound++;
504  h_etaTrack_->Fill(dau1->eta()); //Denominator eta Sta
505  h_ptTrack_->Fill(dau1->pt()); //Denominator pt Sta
506  h_DELTA_ZMuTrackMassReco_dimuonMassGen_->Fill(zMuTrkCand.mass()-dimuonMassGen);
507  }
508 
509  }
510  }
511  } //end loop on Candidate
512 
513  if (!ZMuMuMatchedfound && !ZMuStaMatchedfound && ZMuTrackMatchedfound == 0) noMCmatching++;
514  if (!ZMuMuMatchedfound && ZMuTrackMatchedfound == 1) ZMuTrack_exclusive_1match++;
515  if (!ZMuMuMatchedfound && ZMuTrackMatchedfound > 1) ZMuTrack_exclusive_morematch++;
516  if (!ZMuMuMatchedfound && ZMuTrackMatchedSelected_exclusivefound == 1) ZMuTrackselected_exclusive_1match++;
517  if (!ZMuMuMatchedfound && ZMuTrackMatchedSelected_exclusivefound > 1) ZMuTrackselected_exclusive_morematch++;
518  if (ZMuMuMatchedfound && ZMuTrackMatchedfound == 1) ZMuTrack_ZMuMu_1match++;
519  if (ZMuMuMatchedfound && ZMuTrackMatchedfound == 2) ZMuTrack_ZMuMu_2match++;
520  if (ZMuMuMatchedfound && ZMuTrackMatchedfound > 2) ZMuTrack_ZMuMu_morematch++;
521 
522 }
523 
524 bool ZMuMuEfficiency::check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
525 {
526  int partId0 = dauGen0->pdgId();
527  int partId1 = dauGen1->pdgId();
528  int partId2 = dauGen2->pdgId();
529  bool muplusFound=false;
530  bool muminusFound=false;
531  bool ZFound=false;
532  if (partId0==13 || partId1==13 || partId2==13) muminusFound=true;
533  if (partId0==-13 || partId1==-13 || partId2==-13) muplusFound=true;
534  if (partId0==23 || partId1==23 || partId2==23) ZFound=true;
535  return muplusFound*muminusFound*ZFound;
536 }
537 
538 float ZMuMuEfficiency::getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
539 {
540  int partId0 = dauGen0->pdgId();
541  int partId1 = dauGen1->pdgId();
542  int partId2 = dauGen2->pdgId();
543  float ptpart=0.;
544  if (partId0 == ipart) {
545  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
546  const Candidate * dauMuGen = dauGen0->daughter(k);
547  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
548  ptpart = dauMuGen->pt();
549  }
550  }
551  }
552  if (partId1 == ipart) {
553  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
554  const Candidate * dauMuGen = dauGen1->daughter(k);
555  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
556  ptpart = dauMuGen->pt();
557  }
558  }
559  }
560  if (partId2 == ipart) {
561  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
562  const Candidate * dauMuGen = dauGen2->daughter(k);
563  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
564  ptpart = dauMuGen->pt();
565  }
566  }
567  }
568  return ptpart;
569 }
570 
571 float ZMuMuEfficiency::getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
572 {
573  int partId0 = dauGen0->pdgId();
574  int partId1 = dauGen1->pdgId();
575  int partId2 = dauGen2->pdgId();
576  float etapart=0.;
577  if (partId0 == ipart) {
578  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
579  const Candidate * dauMuGen = dauGen0->daughter(k);
580  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
581  etapart = dauMuGen->eta();
582  }
583  }
584  }
585  if (partId1 == ipart) {
586  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
587  const Candidate * dauMuGen = dauGen1->daughter(k);
588  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
589  etapart = dauMuGen->eta();
590  }
591  }
592  }
593  if (partId2 == ipart) {
594  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
595  const Candidate * dauMuGen = dauGen2->daughter(k);
596  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
597  etapart = dauMuGen->eta();
598  }
599  }
600  }
601  return etapart;
602 }
603 
604 float ZMuMuEfficiency::getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
605 {
606  int partId0 = dauGen0->pdgId();
607  int partId1 = dauGen1->pdgId();
608  int partId2 = dauGen2->pdgId();
609  float phipart=0.;
610  if (partId0 == ipart) {
611  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
612  const Candidate * dauMuGen = dauGen0->daughter(k);
613  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
614  phipart = dauMuGen->phi();
615  }
616  }
617  }
618  if (partId1 == ipart) {
619  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
620  const Candidate * dauMuGen = dauGen1->daughter(k);
621  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
622  phipart = dauMuGen->phi();
623  }
624  }
625  }
626  if (partId2 == ipart) {
627  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
628  const Candidate * dauMuGen = dauGen2->daughter(k);
629  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
630  phipart = dauMuGen->phi();
631  }
632  }
633  }
634  return phipart;
635 }
636 
637 Particle::LorentzVector ZMuMuEfficiency::getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
638 {
639  int partId0 = dauGen0->pdgId();
640  int partId1 = dauGen1->pdgId();
641  int partId2 = dauGen2->pdgId();
642  Particle::LorentzVector p4part(0.,0.,0.,0.);
643  if (partId0 == ipart) {
644  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
645  const Candidate * dauMuGen = dauGen0->daughter(k);
646  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
647  p4part = dauMuGen->p4();
648  }
649  }
650  }
651  if (partId1 == ipart) {
652  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
653  const Candidate * dauMuGen = dauGen1->daughter(k);
654  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
655  p4part = dauMuGen->p4();
656  }
657  }
658  }
659  if (partId2 == ipart) {
660  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
661  const Candidate * dauMuGen = dauGen2->daughter(k);
662  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
663  p4part = dauMuGen->p4();
664  }
665  }
666  }
667  return p4part;
668 }
669 
670 
671 
673  // double efficiencySTA =(double)numberOfOverlappedStandAlone_/(double)numberOfMatchedZMuTrack_;
674  // double errorEff_STA = sqrt( efficiencySTA*(1 - efficiencySTA)/(double)numberOfMatchedZMuTrack_);
675 
677  double myErrTrackEff = sqrt(myTrackEff*(1-myTrackEff)/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuSta_));
678 
680  double myErrStaEff = sqrt(myTrackEff*(1-myTrackEff)/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuTrack_exclusive));
681 
682  // double efficiencyTRACK =(double)numberOfOverlappedTracks_/(double)numberOfMatchedZMuSta_;
683  // double errorEff_TRACK = sqrt( efficiencyTRACK*(1 - efficiencyTRACK)/(double)numberOfMatchedZMuSta_);
684 
685  cout << "------------------------------------ Counters for MC acceptance --------------------------------" << endl;
686  cout << "totalNumberOfevents = " << totalNumberOfevents << endl;
687  cout << "totalNumberOfZfound = " << totalNumberOfZfound << endl;
688  cout << "totalNumberOfZpassed = " << totalNumberOfZPassed << endl;
689  cout << "n. of events zMuMu found (gen level selected)" << n_zMuMufound_genZsele << endl;
690  cout << "n. of events zMuSta found (gen level selected)" << n_zMuStafound_genZsele << endl;
691  cout << "n. of events zMuTrk found (gen level selected)" << n_zMuTrkfound_genZsele << endl;
692 
693  cout << "---------------------------- Counter for MC truth efficiency calculation--------------------- " << endl;
694 
695  cout << "number of events with ZMuMu found = " << numberOfEventsWithZMuMufound << endl;
696  cout << "number of events with ZMuSta found = " << numberOfEventsWithZMuStafound << endl;
697  cout << "-------------------------------------------------------------------------------------- " << endl;
698 
699  cout << "number of events without MC maching = " << noMCmatching << endl;
700  cout << "number of ZMuTrack exclsive 1 match = " << ZMuTrack_exclusive_1match << endl;
701  cout << "number of ZMuTrack exclsive more match = " << ZMuTrack_exclusive_morematch << endl;
702  cout << "number of ZMuTrack selected exclusive 1 match = " << ZMuTrackselected_exclusive_1match << endl;
703  cout << "number of ZMuTrack selected exclusive more match = " << ZMuTrackselected_exclusive_morematch << endl;
704  cout << "number of ZMuTrack ZMuMu 1 match = " << ZMuTrack_ZMuMu_1match << endl;
705  cout << "number of ZMuTrack ZMuMu 2 match = " << ZMuTrack_ZMuMu_2match << endl;
706  cout << "number of ZMuTrack ZMuMu more match = " << ZMuTrack_ZMuMu_morematch << endl;
707  cout << "numberOfMatchedZMuMu = " << numberOfMatchedZMuMu_ << endl;
708  cout << "numberOfMatchedSelectdZMuMu = " << numberOfMatchedSelectedZMuMu_ << endl;
709  cout << "numberOfMatchedZMuSta = " << numberOfMatchedZMuSta_ << endl;
710  cout << "numberOfMatchedSelectedZMuSta = " << numberOfMatchedSelectedZMuSta_ << endl;
711  cout << "numberOfMatchedZMuTrack_matchedZMuMu = " << numberOfMatchedZMuTrack_matchedZMuMu << endl;
712  cout << "numberOfMatchedZMuTrack_matchedSelectedZMuMu = " << numberOfMatchedZMuTrack_matchedSelectedZMuMu << endl;
713  cout << "numberOfMatchedZMuTrack exclusive = " << numberOfMatchedZMuTrack_exclusive << endl;
714  cout << "numberOfMatchedSelectedZMuTrack exclusive = " << numberOfMatchedSelectedZMuTrack_exclusive << endl;
715  cout << " ----------------------------- Efficiency --------------------------------- " << endl;
716  cout << "Efficiency StandAlone = " << myStaEff << " +/- " << myErrStaEff << endl;
717  cout << "Efficiency Track = " << myTrackEff << " +/- " << myErrTrackEff << endl;
718 }
719 
721 
723 
int i
Definition: DBlmapReader.cc:9
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
TH1D * h_ptMuonOverlappedToTrack_
TH1D * h_dimuonPtGenPassed_
ZMuMuEfficiency(const edm::ParameterSet &pset)
edm::InputTag muonIso_
const std::vector< reco::PFCandidatePtr > & tracks_
virtual float mass() const =0
mass
bool check_ifZmumu(const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
int numberOfMatchedSelectedZMuSta_
virtual float eta() const =0
momentum pseudorapidity
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
TH1D * h_etaMuonOverlappedToStandAlone_
virtual int status() const =0
status word
edm::InputTag tracks_
edm::InputTag zMuTrackMatchMap_
int ZMuTrackselected_exclusive_1match
edm::InputTag genParticles_
int numberOfMatchedZMuTrack_exclusive
virtual float phi() const =0
momentum azimuthal angle
TH1D * h_dimuonEtaGenPassed_
int numberOfMatchedZMuTrack_matchedSelectedZMuMu
reco::CandidateBaseRef standAloneMuonCandRef_
int numberOfOverlappedStandAlone_
TH1D * h_DELTA_ZMuMuMassReco_dimuonMassGen_
int numberOfMatchedZMuTrack_matchedZMuMu
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
virtual size_type numberOfDaughters() const =0
number of daughters
virtual float pt() const =0
transverse momentum
edm::ValueMap< float > IsolationCollection
OverlapChecker overlap_
float getParticlePt(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
edm::InputTag muonMatchMap_
edm::InputTag zMuMu_
T sqrt(T t)
Definition: SSEVec.h:48
edm::InputTag zMuMuMatchMap_
virtual void endJob() override
edm::InputTag standAlone_
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
int ZMuTrackselected_exclusive_morematch
int numberOfMatchedSelectedZMuTrack_exclusive
bool isMatched(TrackingRecHit const &hit)
float getParticleEta(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
T * make(const Args &...args) const
make new ROOT object
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
reco::CandidateBaseRef trackMuonCandRef_
edm::InputTag trackIso_
float getParticlePhi(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
edm::InputTag zMuStandAlone_
int k[5][pyjets_maxn]
TH1D * h_DELTA_ZMuStaMassReco_dimuonMassGen_
virtual int pdgId() const =0
PDG identifier.
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
virtual void analyze(const edm::Event &event, const edm::EventSetup &setup) override
tuple tracks
Definition: testEve_cfg.py:39
int numberOfEventsWithZMuMufound
int numberOfMatchedZMuTrack_notOverlapped
unsigned int nbinsPt_
Particle::LorentzVector getParticleP4(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
int numberOfMatchedSelectedZMuMu_
TH1D * h_dimuonMassGenPassed_
TH1D * h_ptMuonOverlappedToStandAlone_
edm::InputTag muons_
tuple muons
Definition: patZpeak.py:38
edm::InputTag zMuStandAloneMatchMap_
TH1D * h_DELTA_ZMuTrackMassReco_dimuonMassGen_
edm::InputTag zMuTrack_
tuple cout
Definition: gather_cfg.py:121
int numberOfEventsWithZMuStafound
TH1D * h_etaMuonOverlappedToTrack_
int ZMuTrack_exclusive_morematch
unsigned int nbinsEta_
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:27
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
edm::InputTag standAloneIso_
tuple zMuMu
zMuMu vector of PSet is common to all categories except zMuTrk category
reco::CandidateBaseRef globalMuonCandRef_
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
virtual bool isGlobalMuon() const =0
virtual const CandidateBaseRef & masterClone() const =0