CMS 3D CMS Logo

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 
27 #include "TH1.h"
28 #include <vector>
29 
30 using namespace edm;
31 using namespace std;
32 using namespace reco;
33 
35 
37 public:
39 private:
40  void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
41  bool check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
42  float getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
43  float getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
44  float getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
45  Particle::LorentzVector getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2);
46  void endJob() override;
47 
62 
63  double zMassMin_, zMassMax_, ptmin_, etamax_, isomax_;
64  unsigned int nbinsPt_, nbinsEta_;
65  reco::CandidateBaseRef globalMuonCandRef_, trackMuonCandRef_, standAloneMuonCandRef_;
67 
68  //histograms for measuring tracker efficiency
69  TH1D *h_etaStandAlone_, *h_etaMuonOverlappedToStandAlone_;
70  TH1D *h_ptStandAlone_, *h_ptMuonOverlappedToStandAlone_;
71 
72  //histograms for measuring standalone efficiency
73  TH1D *h_etaTrack_, *h_etaMuonOverlappedToTrack_;
74  TH1D *h_ptTrack_, *h_ptMuonOverlappedToTrack_;
75 
76  //histograms for MC acceptance
77  TH1D *h_nZMCfound_;
78  TH1D *h_ZetaGen_, *h_ZptGen_, *h_ZmassGen_;
79  TH1D *h_muetaGen_, *h_muptGen_, *h_muIsoGen_;
80  TH1D *h_dimuonPtGen_, *h_dimuonMassGen_, *h_dimuonEtaGen_;
81  TH1D *h_ZetaGenPassed_, *h_ZptGenPassed_, *h_ZmassGenPassed_;
82  TH1D *h_muetaGenPassed_, *h_muptGenPassed_, *h_muIsoGenPassed_;
83  TH1D *h_dimuonPtGenPassed_, *h_dimuonMassGenPassed_, *h_dimuonEtaGenPassed_;
84  //histograms for invarian mass resolution
85  TH1D *h_DELTA_ZMuMuMassReco_dimuonMassGen_, *h_DELTA_ZMuStaMassReco_dimuonMassGen_, *h_DELTA_ZMuTrackMassReco_dimuonMassGen_;
86 
87  int numberOfEventsWithZMuMufound, numberOfEventsWithZMuStafound;
88  int numberOfMatchedZMuSta_,numberOfMatchedSelectedZMuSta_;
89  int numberOfMatchedZMuMu_, numberOfMatchedSelectedZMuMu_;
90  int numberOfOverlappedStandAlone_, numberOfOverlappedTracks_, numberOfMatchedZMuTrack_notOverlapped;
91  int numberOfMatchedZMuTrack_exclusive ,numberOfMatchedSelectedZMuTrack_exclusive;
92  int numberOfMatchedZMuTrack_matchedZMuMu, numberOfMatchedZMuTrack_matchedSelectedZMuMu ;
93  int totalNumberOfevents, totalNumberOfZfound, totalNumberOfZPassed;
94  int noMCmatching, ZMuTrack_exclusive_1match, ZMuTrack_exclusive_morematch;
95  int ZMuTrackselected_exclusive_1match, ZMuTrackselected_exclusive_morematch;
96  int ZMuTrack_ZMuMu_1match, ZMuTrack_ZMuMu_2match, ZMuTrack_ZMuMu_morematch;
97 
98  int n_zMuMufound_genZsele, n_zMuStafound_genZsele, n_zMuTrkfound_genZsele;
99 };
100 
104 #include <iostream>
105 #include <iterator>
106 #include <cmath>
107 
108 
110  zMuMuToken_(consumes<CandidateView>(pset.getParameter<InputTag>("zMuMu"))),
111  zMuMuMatchMapToken_(mayConsume<GenParticleMatch>(pset.getParameter<InputTag>("zMuMuMatchMap"))),
112  zMuTrackToken_(consumes<CandidateView>(pset.getParameter<InputTag>("zMuTrack"))),
113  zMuTrackMatchMapToken_(mayConsume<GenParticleMatch>(pset.getParameter<InputTag>("zMuTrackMatchMap"))),
114  zMuStandAloneToken_(consumes<CandidateView>(pset.getParameter<InputTag>("zMuStandAlone"))),
115  zMuStandAloneMatchMapToken_(mayConsume<GenParticleMatch>(pset.getParameter<InputTag>("zMuStandAloneMatchMap"))),
116  muonsToken_(consumes<CandidateView>(pset.getParameter<InputTag>("muons"))),
117  muonMatchMapToken_(mayConsume<GenParticleMatch>(pset.getParameter<InputTag>("muonMatchMap"))),
118  muonIsoToken_(mayConsume<IsolationCollection>(pset.getParameter<InputTag>("muonIso"))),
119  tracksToken_(consumes<CandidateView>(pset.getParameter<InputTag>("tracks"))),
120  trackIsoToken_(mayConsume<IsolationCollection>(pset.getParameter<InputTag>("trackIso"))),
121  standAloneToken_(consumes<CandidateView>(pset.getParameter<InputTag>("standAlone"))),
122  standAloneIsoToken_(mayConsume<IsolationCollection>(pset.getParameter<InputTag>("standAloneIso"))),
123  genParticlesToken_(consumes<GenParticleCollection>(pset.getParameter<InputTag>( "genParticles"))),
124 
125  zMassMin_(pset.getUntrackedParameter<double>("zMassMin")),
126  zMassMax_(pset.getUntrackedParameter<double>("zMassMax")),
127  ptmin_(pset.getUntrackedParameter<double>("ptmin")),
128  etamax_(pset.getUntrackedParameter<double>("etamax")),
129  isomax_(pset.getUntrackedParameter<double>("isomax")),
130  nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
131  nbinsEta_(pset.getUntrackedParameter<unsigned int>("nbinsEta")) {
133  TFileDirectory trackEffDir = fs->mkdir("TrackEfficiency");
134 
135  // tracker efficiency distributions
136  h_etaStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonEta",
137  "StandAlone #eta for Z -> #mu + standalone",
139  h_etaMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAloneEta",
140  "Global muon overlapped to standAlone #eta for Z -> #mu + sa",
142  h_ptStandAlone_ = trackEffDir.make<TH1D>("StandAloneMuonPt",
143  "StandAlone p_{t} for Z -> #mu + standalone",
144  nbinsPt_, ptmin_, 100);
145  h_ptMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAlonePt",
146  "Global muon overlapped to standAlone p_{t} for Z -> #mu + sa",
147  nbinsPt_, ptmin_, 100);
148 
149 
150  // StandAlone efficiency distributions
151  TFileDirectory standaloneEffDir = fs->mkdir("StandaloneEfficiency");
152  h_etaTrack_ = standaloneEffDir.make<TH1D>("TrackMuonEta",
153  "Track #eta for Z -> #mu + track",
155  h_etaMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackEta",
156  "Global muon overlapped to track #eta for Z -> #mu + tk",
158  h_ptTrack_ = standaloneEffDir.make<TH1D>("TrackMuonPt",
159  "Track p_{t} for Z -> #mu + track",
160  nbinsPt_, ptmin_, 100);
161  h_ptMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackPt",
162  "Global muon overlapped to track p_{t} for Z -> #mu + tk",
163  nbinsPt_, ptmin_, 100);
164 
165 
166  // inv. mass resolution studies
167  TFileDirectory invMassResolutionDir = fs->mkdir("invriantMassResolution");
168  h_DELTA_ZMuMuMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuMu_invMassResolution","zMuMu invariant Mass Resolution",50,-25,25);
169  h_DELTA_ZMuStaMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuSta_invMassResolution","zMuSta invariant Mass Resolution",50,-25,25);
170  h_DELTA_ZMuTrackMassReco_dimuonMassGen_ = invMassResolutionDir.make<TH1D>("zMuTrack_invMassResolution","zMuTrack invariant Mass Resolution",50,-25,25);
171 
172 
173  // generator level histograms
174  TFileDirectory genParticleDir = fs->mkdir("genParticle");
175  h_nZMCfound_ = genParticleDir.make<TH1D>("NumberOfgeneratedZeta","n. of generated Z per event",4,-.5,3.5);
176  h_ZetaGen_ = genParticleDir.make<TH1D>("generatedZeta","#eta of generated Z",100,-5.,5.);
177  h_ZptGen_ = genParticleDir.make<TH1D>("generatedZpt","pt of generated Z",100,0.,200.);
178  h_ZmassGen_ = genParticleDir.make<TH1D>("generatedZmass","mass of generated Z",100,0.,200.);
179  h_muetaGen_ = genParticleDir.make<TH1D>("generatedMuonEta","#eta of generated muons from Z decay",100,-5.,5.);
180  h_muptGen_ = genParticleDir.make<TH1D>("generatedMuonpt","pt of generated muons from Z decay",100,0.,200.);
181  h_dimuonEtaGen_ = genParticleDir.make<TH1D>("generatedDimuonEta","#eta of generated dimuon",100,-5.,5.);
182  h_dimuonPtGen_ = genParticleDir.make<TH1D>("generatedDimuonPt","pt of generated dimuon",100,0.,200.);
183  h_dimuonMassGen_ = genParticleDir.make<TH1D>("generatedDimuonMass","mass of generated dimuon",100,0.,200.);
184  h_ZetaGenPassed_ = genParticleDir.make<TH1D>("generatedZeta_passed","#eta of generated Z after cuts",100,-5.,5.);
185  h_ZptGenPassed_ = genParticleDir.make<TH1D>("generatedZpt_passed","pt of generated Z after cuts",100,0.,200.);
186  h_ZmassGenPassed_ = genParticleDir.make<TH1D>("generatedZmass_passed","mass of generated Z after cuts",100,0.,200.);
187  h_muetaGenPassed_ = genParticleDir.make<TH1D>("generatedMuonEta_passed","#eta of generated muons from Z decay after cuts",100,-5.,5.);
188  h_muptGenPassed_ = genParticleDir.make<TH1D>("generatedMuonpt_passed","pt of generated muons from Z decay after cuts",100,0.,200.);
189  h_dimuonEtaGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonEta_passed","#eta of generated dimuon after cuts",100,-5.,5.);
190  h_dimuonPtGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonPt_passed","pt of generated dimuon after cuts",100,0.,200.);
191  h_dimuonMassGenPassed_ = genParticleDir.make<TH1D>("generatedDimuonMass_passed","mass of generated dimuon after cuts",100,0.,200.);
192  // to insert isolation histograms ..............
193 
207  noMCmatching = 0;
215 
219 
220  // generator counters
224 
225 }
226 
228  Handle<CandidateView> zMuMu;
229  Handle<GenParticleMatch> zMuMuMatchMap; //Map of Z made by Mu global + Mu global
230  Handle<CandidateView> zMuTrack;
231  Handle<GenParticleMatch> zMuTrackMatchMap; //Map of Z made by Mu + Track
232  Handle<CandidateView> zMuStandAlone;
233  Handle<GenParticleMatch> zMuStandAloneMatchMap; //Map of Z made by Mu + StandAlone
234  Handle<CandidateView> muons; //Collection of Muons
235  Handle<GenParticleMatch> muonMatchMap;
237  Handle<CandidateView> tracks; //Collection of Tracks
239  Handle<CandidateView> standAlone; //Collection of StandAlone
240  Handle<IsolationCollection> standAloneIso;
241  Handle<GenParticleCollection> genParticles; // Collection of Generatd Particles
242 
243  event.getByToken(zMuMuToken_, zMuMu);
244  event.getByToken(zMuTrackToken_, zMuTrack);
245  event.getByToken(zMuStandAloneToken_, zMuStandAlone);
246  event.getByToken(muonsToken_, muons);
247  event.getByToken(tracksToken_, tracks);
248  event.getByToken(standAloneToken_, standAlone);
249  event.getByToken(genParticlesToken_, genParticles);
250 
251  cout << "********* zMuMu size : " << zMuMu->size() << endl;
252  cout << "********* zMuStandAlone size : " << zMuStandAlone->size() << endl;
253  cout << "********* zMuTrack size : " << zMuTrack->size() << endl;
254  cout << "********* muons size : " << muons->size()<< endl;
255  cout << "********* standAlone size : " << standAlone->size()<< endl;
256  cout << "********* tracks size : " << tracks->size()<< endl;
257  cout << "********* generated size : " << genParticles->size()<< endl;
258 
259 
260  // generator level distributions
261 
262  int nZMCfound = 0;
264  int ngen = genParticles->size();
265  bool ZMuMuMatchedfound = false;
266  bool ZMuMuMatchedSelectedfound = false;
267  bool ZMuStaMatchedfound = false;
268  //bool ZMuStaMatchedSelectedfound = false;
269  int ZMuTrackMatchedfound = 0;
270  int ZMuTrackMatchedSelected_exclusivefound = 0;
271 
272  double dimuonMassGen = 0;
273 
274  for (int i=0; i<ngen; i++) {
275  const Candidate &genCand = (*genParticles)[i];
276 
277  // if((genCand.pdgId() == 23) && (genCand.status() == 2)) //this is an intermediate Z0
278  // cout << ">>> intermediate Z0 found, with " << genCand.numberOfDaughters() << " daughters" << endl;
279  if((genCand.pdgId() == 23)&&(genCand.status() == 3)) { //this is a Z0
280  if(genCand.numberOfDaughters() == 3) { // possible Z0 decays in mu+ mu-, the 3rd daughter is the same Z0
281  const Candidate * dauGen0 = genCand.daughter(0);
282  const Candidate * dauGen1 = genCand.daughter(1);
283  const Candidate * dauGen2 = genCand.daughter(2);
284  if (check_ifZmumu(dauGen0, dauGen1, dauGen2)) { // Z0 in mu+ mu-
286  nZMCfound++;
287  bool checkpt = false;
288  bool checketa = false;
289  bool checkmass = false;
290  float mupluspt, muminuspt, mupluseta, muminuseta;
291  mupluspt = getParticlePt(-13,dauGen0,dauGen1,dauGen2);
292  muminuspt = getParticlePt(13,dauGen0,dauGen1,dauGen2);
293  mupluseta = getParticleEta(-13,dauGen0,dauGen1,dauGen2);
294  muminuseta = getParticleEta(13,dauGen0,dauGen1,dauGen2);
295  //float muplusphi = getParticlePhi(-13,dauGen0,dauGen1,dauGen2);
296  //float muminusphi = getParticlePhi(13,dauGen0,dauGen1,dauGen2);
297  Particle::LorentzVector pZ(0, 0, 0, 0);
298  Particle::LorentzVector muplusp4 = getParticleP4(-13,dauGen0,dauGen1,dauGen2);
299  Particle::LorentzVector muminusp4 = getParticleP4(13,dauGen0,dauGen1,dauGen2);
300  pZ = muplusp4 + muminusp4;
301  double dimuon_pt = sqrt(pZ.x()*pZ.x()+pZ.y()*pZ.y());
302  double tan_theta_half = tan(atan(dimuon_pt/pZ.z())/2.);
303  double dimuon_eta = 0.;
304  if (tan_theta_half>0) dimuon_eta = -log(tan(tan_theta_half));
305  if (tan_theta_half<=0) dimuon_eta = log(tan(-tan_theta_half));
306 
307  dimuonMassGen = pZ.mass(); // dimuon invariant Mass at Generator Level
308 
309  h_ZmassGen_->Fill(genCand.mass());
310  h_ZetaGen_->Fill(genCand.eta());
311  h_ZptGen_->Fill(genCand.pt());
312  h_dimuonMassGen_->Fill(pZ.mass());
313  h_dimuonEtaGen_->Fill(dimuon_eta);
314  h_dimuonPtGen_->Fill(dimuon_pt);
315  h_muetaGen_->Fill(mupluseta);
316  h_muetaGen_->Fill(muminuseta);
317  h_muptGen_->Fill(mupluspt);
318  h_muptGen_->Fill(muminuspt);
319  // dimuon 4-momentum
320  // h_mDimuonMC->Fill(pZ.mass());
321  // h_ZminusDimuonMassMC->Fill(genCand.mass()-pZ.mass());
322  // h_DeltaPhiMC->Fill(deltaPhi(muplusphi,muminusphi));
323  // if (dauGen2==23) float z_eta = dauGen2->eta();
324  // if (dauGen2==23) float Zpt = dauGen2->pt();
325  // h_DeltaPhivsZPtMC->Fill(DeltaPhi(muplusphi,muminusphi),ZPt);
326 
327  if (mupluspt > ptmin_ && muminuspt > ptmin_) checkpt = true;
328  if (mupluseta < etamax_ && muminuseta < etamax_) checketa = true;
329  if (genCand.mass()>zMassMin_ && genCand.mass()<zMassMax_) checkmass = true;
330  if (checkpt && checketa && checkmass) {
332  h_ZmassGenPassed_->Fill(genCand.mass());
333  h_ZetaGenPassed_->Fill(genCand.eta());
334  h_ZptGenPassed_->Fill(genCand.pt());
335  h_dimuonMassGenPassed_->Fill(pZ.mass());
336  h_dimuonEtaGenPassed_->Fill(dimuon_eta);
337  h_dimuonPtGenPassed_->Fill(dimuon_pt);
338  h_muetaGenPassed_->Fill(mupluseta);
339  h_muetaGenPassed_->Fill(muminuseta);
340  h_muptGenPassed_->Fill(mupluspt);
341  h_muptGenPassed_->Fill(muminuspt);
342 
343  if (!zMuMu->empty() ) {
345  }
346  else if (!zMuStandAlone->empty() ) {
348  }
349  else {
351  }
352 
353  }
354  }
355 
356  }
357  }
358  }
359  h_nZMCfound_->Fill(nZMCfound); // number of Z found in the event at generator level
360 
361  //TRACK efficiency (conto numero di eventi Zmumu global e ZmuSta (ricorda che sono due campioni esclusivi)
362 
363  if (!zMuMu->empty() ) {
365  event.getByToken(zMuMuMatchMapToken_, zMuMuMatchMap);
366  event.getByToken(muonIsoToken_, muonIso);
367  event.getByToken(standAloneIsoToken_, standAloneIso);
368  event.getByToken(muonMatchMapToken_, muonMatchMap);
369  for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates
370  const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate
371  CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i);
372  bool isMatched = false;
373  GenParticleRef zMuMuMatch = (*zMuMuMatchMap)[zMuMuCandRef];
374 
375  if(zMuMuMatch.isNonnull()) { // ZMuMu matched
376  isMatched = true;
378  }
379  CandidateBaseRef dau0 = zMuMuCand.daughter(0)->masterClone();
380  CandidateBaseRef dau1 = zMuMuCand.daughter(1)->masterClone();
381  if (isMatched) ZMuMuMatchedfound = true;
382 
383  // Cuts
384  if((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) &&
385  (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta()) < etamax_) &&
386  (zMuMuCand.mass() > zMassMin_) && (zMuMuCand.mass() < zMassMax_) &&
387  (isMatched)) {
388  //The Z daughters are already matched!
389  const double globalMuonIsolation0 = (*muonIso)[dau0];
390  const double globalMuonIsolation1 = (*muonIso)[dau1];
391  if((globalMuonIsolation0 < isomax_) && (globalMuonIsolation1 < isomax_)) { // ZMuMu matched and selected by cuts
392  ZMuMuMatchedSelectedfound = true;
394  h_etaStandAlone_->Fill(dau0->eta()); // StandAlone found dau0, eta
395  h_etaStandAlone_->Fill(dau1->eta()); // StandAlone found dau1, eta
396  h_etaMuonOverlappedToStandAlone_->Fill(dau0->eta()); // is global muon so dau0 is also found as a track, eta
397  h_etaMuonOverlappedToStandAlone_->Fill(dau1->eta()); // is global muon so dau1 is also found as a track, eta
398  h_ptStandAlone_->Fill(dau0->pt()); // StandAlone found dau0, pt
399  h_ptStandAlone_->Fill(dau1->pt()); // StandAlone found dau1, pt
400  h_ptMuonOverlappedToStandAlone_->Fill(dau0->pt()); // is global muon so dau0 is also found as a track, pt
401  h_ptMuonOverlappedToStandAlone_->Fill(dau1->pt()); // is global muon so dau1 is also found as a track, pt
402 
403  h_etaTrack_->Fill(dau0->eta()); // Track found dau0, eta
404  h_etaTrack_->Fill(dau1->eta()); // Track found dau1, eta
405  h_etaMuonOverlappedToTrack_->Fill(dau0->eta()); // is global muon so dau0 is also found as a StandAlone, eta
406  h_etaMuonOverlappedToTrack_->Fill(dau1->eta()); // is global muon so dau1 is also found as a StandAlone, eta
407  h_ptTrack_->Fill(dau0->pt()); // Track found dau0, pt
408  h_ptTrack_->Fill(dau1->pt()); // Track found dau1, pt
409  h_ptMuonOverlappedToTrack_->Fill(dau0->pt()); // is global muon so dau0 is also found as a StandAlone, pt
410  h_ptMuonOverlappedToTrack_->Fill(dau1->pt()); // is global muon so dau1 is also found as a StandAlone, pt
411 
412  h_DELTA_ZMuMuMassReco_dimuonMassGen_->Fill(zMuMuCand.mass()-dimuonMassGen);
413  // check that the two muons are matched . .per ora è solo un mio controllo
414  for(unsigned int j = 0; j < muons->size() ; ++j) {
415  CandidateBaseRef muCandRef = muons->refAt(j);
416  GenParticleRef muonMatch = (*muonMatchMap)[muCandRef];
417  // if (muonMatch.isNonnull()) cout << "mu match n. " << j << endl;
418  }
419  }
420  }
421  }
422  }
423 
424  if (!zMuStandAlone->empty()) {
426  event.getByToken(zMuStandAloneMatchMapToken_, zMuStandAloneMatchMap);
427  event.getByToken(muonIsoToken_, muonIso);
428  event.getByToken(standAloneIsoToken_, standAloneIso);
429  event.getByToken(muonMatchMapToken_, muonMatchMap);
430  for(unsigned int i = 0; i < zMuStandAlone->size(); ++i) { //loop on candidates
431  const Candidate & zMuStaCand = (*zMuStandAlone)[i]; //the candidate
432  CandidateBaseRef zMuStaCandRef = zMuStandAlone->refAt(i);
433  bool isMatched = false;
434  GenParticleRef zMuStaMatch = (*zMuStandAloneMatchMap)[zMuStaCandRef];
435  if(zMuStaMatch.isNonnull()) { // ZMuSta Macthed
436  isMatched = true;
437  ZMuStaMatchedfound = true;
439  }
440  CandidateBaseRef dau0 = zMuStaCand.daughter(0)->masterClone();
441  CandidateBaseRef dau1 = zMuStaCand.daughter(1)->masterClone();
442 
443  // Cuts
444  if((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) &&
445  (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta()) < etamax_) &&
446  (zMuStaCand.mass() > zMassMin_) && (zMuStaCand.mass() < zMassMax_) &&
447  (isMatched)) {
449  if(dau0->isGlobalMuon()) {
450  standAloneMuonCandRef_ = dau1;
451  globalMuonCandRef_ = dau0;
452  }
453  if(dau1->isGlobalMuon()) {
454  standAloneMuonCandRef_ = dau0;
455  globalMuonCandRef_ = dau1;
456  }
457 
458  const double globalMuonIsolation = (*muonIso)[globalMuonCandRef_];
459  const double standAloneMuonIsolation = (*standAloneIso)[standAloneMuonCandRef_];
460 
461  if((globalMuonIsolation < isomax_) && (standAloneMuonIsolation < isomax_)) { // ZMuSta matched and selected
462  //ZMuStaMatchedSelectedfound = true;
464  h_etaStandAlone_->Fill(standAloneMuonCandRef_->eta()); //Denominator eta for measuring track efficiency
465  h_ptStandAlone_->Fill(standAloneMuonCandRef_->pt()); //Denominator pt for measuring track eff
466  h_DELTA_ZMuStaMassReco_dimuonMassGen_->Fill(zMuStaCand.mass()-dimuonMassGen); // differnce between ZMuSta reco and dimuon mass gen
467 
468  }
469  }
470  }
471  } //end loop on Candidate
472 
473  //STANDALONE efficiency
474 
475  if (!zMuTrack->empty()) {
476  event.getByToken(zMuTrackMatchMapToken_, zMuTrackMatchMap);
477  event.getByToken(muonIsoToken_, muonIso);
478  event.getByToken(trackIsoToken_, trackIso);
479  event.getByToken(muonMatchMapToken_, muonMatchMap);
480  for(unsigned int i = 0; i < zMuTrack->size(); ++i) { //loop on candidates
481  const Candidate & zMuTrkCand = (*zMuTrack)[i]; //the candidate
482  CandidateBaseRef zMuTrkCandRef = zMuTrack->refAt(i);
483  bool isMatched = false;
484  GenParticleRef zMuTrkMatch = (*zMuTrackMatchMap)[zMuTrkCandRef];
485  if(zMuTrkMatch.isNonnull()) {
486  isMatched = true;
487  }
488  CandidateBaseRef dau0 = zMuTrkCand.daughter(0)->masterClone();
489  CandidateBaseRef dau1 = zMuTrkCand.daughter(1)->masterClone();
490 
491  if (isMatched) {
492  ZMuTrackMatchedfound++;
493  if (ZMuMuMatchedfound) numberOfMatchedZMuTrack_matchedZMuMu++;
494  if (ZMuMuMatchedSelectedfound) numberOfMatchedZMuTrack_matchedSelectedZMuMu++;
495  if (!ZMuMuMatchedfound) numberOfMatchedZMuTrack_exclusive++;
496  }
497  // Cuts
498  if ((dau0->pt() > ptmin_) && (dau1->pt() > ptmin_) &&
499  (fabs(dau0->eta()) < etamax_) && (fabs(dau1->eta())< etamax_) &&
500  (zMuTrkCand.mass() > zMassMin_) && (zMuTrkCand.mass() < zMassMax_) &&
501  (isMatched) && !ZMuMuMatchedfound && !ZMuStaMatchedfound ) {
502 
503  // dau0 is always the global muon, dau1 is the track for ZMuTrack collection
504  const double globalMuonIsolation = (*muonIso)[dau0];
505  const double trackMuonIsolation = (*trackIso)[dau1];
506  if((globalMuonIsolation < isomax_) && (trackMuonIsolation < isomax_)) { // ZMuTRack matched - selected without ZMuMu found (exclusive)
508  ZMuTrackMatchedSelected_exclusivefound++;
509  h_etaTrack_->Fill(dau1->eta()); //Denominator eta Sta
510  h_ptTrack_->Fill(dau1->pt()); //Denominator pt Sta
511  h_DELTA_ZMuTrackMassReco_dimuonMassGen_->Fill(zMuTrkCand.mass()-dimuonMassGen);
512  }
513 
514  }
515  }
516  } //end loop on Candidate
517 
518  if (!ZMuMuMatchedfound && !ZMuStaMatchedfound && ZMuTrackMatchedfound == 0) noMCmatching++;
519  if (!ZMuMuMatchedfound && ZMuTrackMatchedfound == 1) ZMuTrack_exclusive_1match++;
520  if (!ZMuMuMatchedfound && ZMuTrackMatchedfound > 1) ZMuTrack_exclusive_morematch++;
521  if (!ZMuMuMatchedfound && ZMuTrackMatchedSelected_exclusivefound == 1) ZMuTrackselected_exclusive_1match++;
522  if (!ZMuMuMatchedfound && ZMuTrackMatchedSelected_exclusivefound > 1) ZMuTrackselected_exclusive_morematch++;
523  if (ZMuMuMatchedfound && ZMuTrackMatchedfound == 1) ZMuTrack_ZMuMu_1match++;
524  if (ZMuMuMatchedfound && ZMuTrackMatchedfound == 2) ZMuTrack_ZMuMu_2match++;
525  if (ZMuMuMatchedfound && ZMuTrackMatchedfound > 2) ZMuTrack_ZMuMu_morematch++;
526 
527 }
528 
529 bool ZMuMuEfficiency::check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
530 {
531  int partId0 = dauGen0->pdgId();
532  int partId1 = dauGen1->pdgId();
533  int partId2 = dauGen2->pdgId();
534  bool muplusFound=false;
535  bool muminusFound=false;
536  bool ZFound=false;
537  if (partId0==13 || partId1==13 || partId2==13) muminusFound=true;
538  if (partId0==-13 || partId1==-13 || partId2==-13) muplusFound=true;
539  if (partId0==23 || partId1==23 || partId2==23) ZFound=true;
540  return ( muplusFound && muminusFound && ZFound );
541 }
542 
543 float ZMuMuEfficiency::getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
544 {
545  int partId0 = dauGen0->pdgId();
546  int partId1 = dauGen1->pdgId();
547  int partId2 = dauGen2->pdgId();
548  float ptpart=0.;
549  if (partId0 == ipart) {
550  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
551  const Candidate * dauMuGen = dauGen0->daughter(k);
552  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
553  ptpart = dauMuGen->pt();
554  }
555  }
556  }
557  if (partId1 == ipart) {
558  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
559  const Candidate * dauMuGen = dauGen1->daughter(k);
560  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
561  ptpart = dauMuGen->pt();
562  }
563  }
564  }
565  if (partId2 == ipart) {
566  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
567  const Candidate * dauMuGen = dauGen2->daughter(k);
568  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
569  ptpart = dauMuGen->pt();
570  }
571  }
572  }
573  return ptpart;
574 }
575 
576 float ZMuMuEfficiency::getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
577 {
578  int partId0 = dauGen0->pdgId();
579  int partId1 = dauGen1->pdgId();
580  int partId2 = dauGen2->pdgId();
581  float etapart=0.;
582  if (partId0 == ipart) {
583  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
584  const Candidate * dauMuGen = dauGen0->daughter(k);
585  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
586  etapart = dauMuGen->eta();
587  }
588  }
589  }
590  if (partId1 == ipart) {
591  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
592  const Candidate * dauMuGen = dauGen1->daughter(k);
593  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
594  etapart = dauMuGen->eta();
595  }
596  }
597  }
598  if (partId2 == ipart) {
599  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
600  const Candidate * dauMuGen = dauGen2->daughter(k);
601  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
602  etapart = dauMuGen->eta();
603  }
604  }
605  }
606  return etapart;
607 }
608 
609 float ZMuMuEfficiency::getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
610 {
611  int partId0 = dauGen0->pdgId();
612  int partId1 = dauGen1->pdgId();
613  int partId2 = dauGen2->pdgId();
614  float phipart=0.;
615  if (partId0 == ipart) {
616  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
617  const Candidate * dauMuGen = dauGen0->daughter(k);
618  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
619  phipart = dauMuGen->phi();
620  }
621  }
622  }
623  if (partId1 == ipart) {
624  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
625  const Candidate * dauMuGen = dauGen1->daughter(k);
626  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
627  phipart = dauMuGen->phi();
628  }
629  }
630  }
631  if (partId2 == ipart) {
632  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
633  const Candidate * dauMuGen = dauGen2->daughter(k);
634  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
635  phipart = dauMuGen->phi();
636  }
637  }
638  }
639  return phipart;
640 }
641 
642 Particle::LorentzVector ZMuMuEfficiency::getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
643 {
644  int partId0 = dauGen0->pdgId();
645  int partId1 = dauGen1->pdgId();
646  int partId2 = dauGen2->pdgId();
647  Particle::LorentzVector p4part(0.,0.,0.,0.);
648  if (partId0 == ipart) {
649  for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
650  const Candidate * dauMuGen = dauGen0->daughter(k);
651  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
652  p4part = dauMuGen->p4();
653  }
654  }
655  }
656  if (partId1 == ipart) {
657  for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
658  const Candidate * dauMuGen = dauGen1->daughter(k);
659  if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
660  p4part = dauMuGen->p4();
661  }
662  }
663  }
664  if (partId2 == ipart) {
665  for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
666  const Candidate * dauMuGen = dauGen2->daughter(k);
667  if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
668  p4part = dauMuGen->p4();
669  }
670  }
671  }
672  return p4part;
673 }
674 
675 
676 
678  // double efficiencySTA =(double)numberOfOverlappedStandAlone_/(double)numberOfMatchedZMuTrack_;
679  // double errorEff_STA = sqrt( efficiencySTA*(1 - efficiencySTA)/(double)numberOfMatchedZMuTrack_);
680 
682  double myErrTrackEff = sqrt(myTrackEff*(1-myTrackEff)/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuSta_));
683 
685  double myErrStaEff = sqrt(myTrackEff*(1-myTrackEff)/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuTrack_exclusive));
686 
687  // double efficiencyTRACK =(double)numberOfOverlappedTracks_/(double)numberOfMatchedZMuSta_;
688  // double errorEff_TRACK = sqrt( efficiencyTRACK*(1 - efficiencyTRACK)/(double)numberOfMatchedZMuSta_);
689 
690  cout << "------------------------------------ Counters for MC acceptance --------------------------------" << endl;
691  cout << "totalNumberOfevents = " << totalNumberOfevents << endl;
692  cout << "totalNumberOfZfound = " << totalNumberOfZfound << endl;
693  cout << "totalNumberOfZpassed = " << totalNumberOfZPassed << endl;
694  cout << "n. of events zMuMu found (gen level selected)" << n_zMuMufound_genZsele << endl;
695  cout << "n. of events zMuSta found (gen level selected)" << n_zMuStafound_genZsele << endl;
696  cout << "n. of events zMuTrk found (gen level selected)" << n_zMuTrkfound_genZsele << endl;
697 
698  cout << "---------------------------- Counter for MC truth efficiency calculation--------------------- " << endl;
699 
700  cout << "number of events with ZMuMu found = " << numberOfEventsWithZMuMufound << endl;
701  cout << "number of events with ZMuSta found = " << numberOfEventsWithZMuStafound << endl;
702  cout << "-------------------------------------------------------------------------------------- " << endl;
703 
704  cout << "number of events without MC maching = " << noMCmatching << endl;
705  cout << "number of ZMuTrack exclsive 1 match = " << ZMuTrack_exclusive_1match << endl;
706  cout << "number of ZMuTrack exclsive more match = " << ZMuTrack_exclusive_morematch << endl;
707  cout << "number of ZMuTrack selected exclusive 1 match = " << ZMuTrackselected_exclusive_1match << endl;
708  cout << "number of ZMuTrack selected exclusive more match = " << ZMuTrackselected_exclusive_morematch << endl;
709  cout << "number of ZMuTrack ZMuMu 1 match = " << ZMuTrack_ZMuMu_1match << endl;
710  cout << "number of ZMuTrack ZMuMu 2 match = " << ZMuTrack_ZMuMu_2match << endl;
711  cout << "number of ZMuTrack ZMuMu more match = " << ZMuTrack_ZMuMu_morematch << endl;
712  cout << "numberOfMatchedZMuMu = " << numberOfMatchedZMuMu_ << endl;
713  cout << "numberOfMatchedSelectdZMuMu = " << numberOfMatchedSelectedZMuMu_ << endl;
714  cout << "numberOfMatchedZMuSta = " << numberOfMatchedZMuSta_ << endl;
715  cout << "numberOfMatchedSelectedZMuSta = " << numberOfMatchedSelectedZMuSta_ << endl;
716  cout << "numberOfMatchedZMuTrack_matchedZMuMu = " << numberOfMatchedZMuTrack_matchedZMuMu << endl;
717  cout << "numberOfMatchedZMuTrack_matchedSelectedZMuMu = " << numberOfMatchedZMuTrack_matchedSelectedZMuMu << endl;
718  cout << "numberOfMatchedZMuTrack exclusive = " << numberOfMatchedZMuTrack_exclusive << endl;
719  cout << "numberOfMatchedSelectedZMuTrack exclusive = " << numberOfMatchedSelectedZMuTrack_exclusive << endl;
720  cout << " ----------------------------- Efficiency --------------------------------- " << endl;
721  cout << "Efficiency StandAlone = " << myStaEff << " +/- " << myErrStaEff << endl;
722  cout << "Efficiency Track = " << myTrackEff << " +/- " << myErrTrackEff << endl;
723 }
724 
726 
728 
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
EDGetTokenT< GenParticleMatch > zMuMuMatchMapToken_
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:251
TH1D * h_ptMuonOverlappedToTrack_
TH1D * h_dimuonPtGenPassed_
ZMuMuEfficiency(const edm::ParameterSet &pset)
bool check_ifZmumu(const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
EDGetTokenT< CandidateView > zMuStandAloneToken_
int numberOfMatchedSelectedZMuSta_
TH1D * h_etaMuonOverlappedToStandAlone_
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
int ZMuTrackselected_exclusive_1match
int numberOfMatchedZMuTrack_exclusive
def setup(process, global_tag, zero_tesla=False)
Definition: GeneralSetup.py:2
size_type size() const
TH1D * h_dimuonEtaGenPassed_
EDGetTokenT< CandidateView > zMuTrackToken_
int numberOfMatchedZMuTrack_matchedSelectedZMuMu
reco::CandidateBaseRef standAloneMuonCandRef_
int numberOfOverlappedStandAlone_
example_stream void analyze(const edm::Event &, const edm::EventSetup &) override
TH1D * h_DELTA_ZMuMuMassReco_dimuonMassGen_
int numberOfMatchedZMuTrack_matchedZMuMu
RefToBase< value_type > refAt(size_type i) const
EDGetTokenT< CandidateView > tracksToken_
virtual int status() const =0
status word
EDGetTokenT< GenParticleCollection > genParticlesToken_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
OverlapChecker overlap_
float getParticlePt(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
virtual int pdgId() const =0
PDG identifier.
T sqrt(T t)
Definition: SSEVec.h:18
void endJob() override
bool empty() const
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
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
reco::CandidateBaseRef trackMuonCandRef_
EDGetTokenT< CandidateView > zMuMuToken_
float getParticlePhi(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
int k[5][pyjets_maxn]
TH1D * h_DELTA_ZMuStaMassReco_dimuonMassGen_
virtual const CandidateBaseRef & masterClone() const =0
virtual bool isGlobalMuon() const =0
virtual double eta() const =0
momentum pseudorapidity
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
void analyze(const edm::Event &event, const edm::EventSetup &setup) override
virtual double pt() const =0
transverse momentum
int numberOfEventsWithZMuMufound
virtual double mass() const =0
mass
int numberOfMatchedZMuTrack_notOverlapped
unsigned int nbinsPt_
Particle::LorentzVector getParticleP4(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
int numberOfMatchedSelectedZMuMu_
ValueMap< float > IsolationCollection
TH1D * h_dimuonMassGenPassed_
TH1D * h_ptMuonOverlappedToStandAlone_
EDGetTokenT< CandidateView > standAloneToken_
fixed size matrix
HLT enums.
TH1D * h_DELTA_ZMuTrackMassReco_dimuonMassGen_
EDGetTokenT< IsolationCollection > standAloneIsoToken_
EDGetTokenT< GenParticleMatch > zMuTrackMatchMapToken_
int numberOfEventsWithZMuStafound
TH1D * h_etaMuonOverlappedToTrack_
int ZMuTrack_exclusive_morematch
unsigned int nbinsEta_
EDGetTokenT< GenParticleMatch > muonMatchMapToken_
virtual size_type numberOfDaughters() const =0
number of daughters
EDGetTokenT< IsolationCollection > trackIsoToken_
virtual double phi() const =0
momentum azimuthal angle
EDGetTokenT< CandidateView > muonsToken_
EDGetTokenT< IsolationCollection > muonIsoToken_
Definition: event.py:1
math::PtEtaPhiELorentzVectorF LorentzVector
EDGetTokenT< GenParticleMatch > zMuStandAloneMatchMapToken_
reco::CandidateBaseRef globalMuonCandRef_