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 
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  virtual 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  virtual 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 
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->size() > 0 ) {
345  }
346  else if (zMuStandAlone->size() > 0 ) {
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->size() > 0 ) {
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->size() > 0) {
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->size() > 0) {
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
int i
Definition: DBlmapReader.cc:9
EDGetTokenT< GenParticleMatch > zMuMuMatchMapToken_
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)
virtual float mass() const =0
mass
bool check_ifZmumu(const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
EDGetTokenT< CandidateView > zMuStandAloneToken_
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
int ZMuTrackselected_exclusive_1match
int numberOfMatchedZMuTrack_exclusive
virtual float phi() const =0
momentum azimuthal angle
TH1D * h_dimuonEtaGenPassed_
EDGetTokenT< CandidateView > zMuTrackToken_
int numberOfMatchedZMuTrack_matchedSelectedZMuMu
reco::CandidateBaseRef standAloneMuonCandRef_
int numberOfOverlappedStandAlone_
TH1D * h_DELTA_ZMuMuMassReco_dimuonMassGen_
int numberOfMatchedZMuTrack_matchedZMuMu
EDGetTokenT< CandidateView > tracksToken_
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
EDGetTokenT< GenParticleCollection > genParticlesToken_
edm::ValueMap< float > IsolationCollection
OverlapChecker overlap_
float getParticlePt(const int ipart, const Candidate *dauGen0, const Candidate *dauGen1, const Candidate *dauGen2)
T sqrt(T t)
Definition: SSEVec.h:48
virtual void endJob() override
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_
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 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_
EDGetTokenT< CandidateView > standAloneToken_
tuple muons
Definition: patZpeak.py:38
TH1D * h_DELTA_ZMuTrackMassReco_dimuonMassGen_
EDGetTokenT< IsolationCollection > standAloneIsoToken_
EDGetTokenT< GenParticleMatch > zMuTrackMatchMapToken_
tuple cout
Definition: gather_cfg.py:121
int numberOfEventsWithZMuStafound
TH1D * h_etaMuonOverlappedToTrack_
int ZMuTrack_exclusive_morematch
unsigned int nbinsEta_
EDGetTokenT< GenParticleMatch > muonMatchMapToken_
EDGetTokenT< IsolationCollection > trackIsoToken_
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
EDGetTokenT< CandidateView > muonsToken_
EDGetTokenT< IsolationCollection > muonIsoToken_
tuple zMuMu
zMuMu vector of PSet is common to all categories except zMuTrk category
math::PtEtaPhiELorentzVectorF LorentzVector
EDGetTokenT< GenParticleMatch > zMuStandAloneMatchMapToken_
reco::CandidateBaseRef globalMuonCandRef_
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
virtual bool isGlobalMuon() const =0
virtual const CandidateBaseRef & masterClone() const =0