PageRenderTime 113ms CodeModel.GetById 14ms app.highlight 90ms RepoModel.GetById 2ms app.codeStats 0ms

/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMuEfficiency.cc

https://github.com/aivanov-cern/cmssw
C++ | 723 lines | 601 code | 76 blank | 46 comment | 171 complexity | e167282e0fc44917350a41f57757cd1f MD5 | raw 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
 11#include "DataFormats/Common/interface/AssociationVector.h"
 12#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 13#include "DataFormats/Candidate/interface/CandMatchMap.h" 
 14#include "FWCore/Framework/interface/EDAnalyzer.h"
 15#include "DataFormats/Candidate/interface/Candidate.h"
 16#include "DataFormats/Candidate/interface/CandidateFwd.h"
 17#include "FWCore/ParameterSet/interface/ParameterSet.h"
 18#include "FWCore/Framework/interface/Event.h"
 19#include "FWCore/Framework/interface/EventSetup.h"
 20#include "FWCore/Utilities/interface/InputTag.h"
 21#include "DataFormats/Candidate/interface/OverlapChecker.h"
 22#include "TH1.h"
 23#include <vector>
 24using namespace edm;
 25using namespace std;
 26using namespace reco;
 27
 28class ZMuMuEfficiency : public edm::EDAnalyzer {
 29public:
 30  ZMuMuEfficiency(const edm::ParameterSet& pset);
 31private:
 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
 40  edm::InputTag zMuMu_, zMuMuMatchMap_; 
 41  edm::InputTag zMuTrack_, zMuTrackMatchMap_; 
 42  edm::InputTag zMuStandAlone_, zMuStandAloneMatchMap_;
 43  edm::InputTag muons_, muonMatchMap_, muonIso_;
 44  edm::InputTag tracks_, trackIso_;
 45  edm::InputTag standAlone_, standAloneIso_;
 46  edm::InputTag genParticles_;
 47
 48  double zMassMin_, zMassMax_, ptmin_, etamax_, isomax_;
 49  unsigned int nbinsPt_, nbinsEta_;
 50  reco::CandidateBaseRef globalMuonCandRef_, trackMuonCandRef_, standAloneMuonCandRef_;
 51  OverlapChecker overlap_;
 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
 86#include "FWCore/ServiceRegistry/interface/Service.h"
 87#include "CommonTools/UtilAlgos/interface/TFileService.h"
 88#include "DataFormats/Common/interface/Handle.h"
 89#include "DataFormats/Candidate/interface/Particle.h"
 90#include "DataFormats/Candidate/interface/Candidate.h"
 91#include "DataFormats/Common/interface/ValueMap.h"
 92#include "DataFormats/Candidate/interface/CandAssociation.h"
 93#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 94#include <iostream>
 95#include <iterator>
 96#include <cmath>
 97using namespace std;
 98using namespace reco;
 99using namespace edm;
100
101
102typedef edm::ValueMap<float> IsolationCollection;
103
104ZMuMuEfficiency::ZMuMuEfficiency(const ParameterSet& pset) : 
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")) {
127  Service<TFileService> fs;
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", 
133					    nbinsEta_, -etamax_, etamax_);
134  h_etaMuonOverlappedToStandAlone_ = trackEffDir.make<TH1D>("MuonOverlappedToStandAloneEta", 
135							    "Global muon overlapped to standAlone #eta for Z -> #mu + sa", 
136							    nbinsEta_, -etamax_, etamax_);
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", 
149					    nbinsEta_, -etamax_, etamax_);
150  h_etaMuonOverlappedToTrack_ = standaloneEffDir.make<TH1D>("MuonOverlappedToTrackEta", 
151							    "Global muon overlapped to track #eta for Z -> #mu + tk", 
152							    nbinsEta_, -etamax_, etamax_);
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
189  numberOfEventsWithZMuMufound = 0;
190  numberOfEventsWithZMuStafound = 0;
191  numberOfMatchedZMuMu_ = 0;
192  numberOfMatchedSelectedZMuMu_ = 0;
193  numberOfMatchedZMuSta_ = 0;
194  numberOfMatchedSelectedZMuSta_ = 0;
195  numberOfMatchedZMuTrack_matchedZMuMu = 0;
196  numberOfMatchedZMuTrack_matchedSelectedZMuMu = 0;
197  numberOfMatchedZMuTrack_exclusive = 0;
198  numberOfMatchedSelectedZMuTrack_exclusive = 0;
199  numberOfOverlappedStandAlone_ = 0;
200  numberOfOverlappedTracks_ = 0;
201  numberOfMatchedZMuTrack_notOverlapped = 0;
202  noMCmatching = 0;
203  ZMuTrack_exclusive_1match = 0;
204  ZMuTrack_exclusive_morematch = 0;
205  ZMuTrackselected_exclusive_1match = 0;
206  ZMuTrackselected_exclusive_morematch = 0;
207  ZMuTrack_ZMuMu_1match = 0;
208  ZMuTrack_ZMuMu_2match = 0;
209  ZMuTrack_ZMuMu_morematch = 0;
210
211  n_zMuMufound_genZsele = 0;
212  n_zMuStafound_genZsele = 0;
213  n_zMuTrkfound_genZsele = 0;
214
215  // generator counters
216  totalNumberOfevents = 0;
217  totalNumberOfZfound = 0;
218  totalNumberOfZPassed = 0;
219  
220}
221
222void ZMuMuEfficiency::analyze(const Event& event, const EventSetup& setup) {
223  Handle<CandidateView> zMuMu;  
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; 
231  Handle<IsolationCollection> muonIso; 
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;
258  totalNumberOfevents++;
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-
280	  totalNumberOfZfound++;
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) {
326	    totalNumberOfZPassed++;
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 ) {
339	      n_zMuMufound_genZsele++; 
340	    }
341	    else if (zMuStandAlone->size() > 0 ) {
342		n_zMuStafound_genZsele++;
343	    }
344	    else {
345	      n_zMuTrkfound_genZsele++;	      
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 ) {
359    numberOfEventsWithZMuMufound++;
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;
372	numberOfMatchedZMuMu_++;
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;
388	  numberOfMatchedSelectedZMuMu_++;
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) {
420    numberOfEventsWithZMuStafound++;
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;
433	numberOfMatchedZMuSta_++;
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)) {
443	CandidateBaseRef standAloneMuonCandRef_, globalMuonCandRef_;
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;
458	  numberOfMatchedSelectedZMuSta_++;
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)
502	  numberOfMatchedSelectedZMuTrack_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
524bool 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 
538float 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 
571float 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
604float 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
637Particle::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
672void ZMuMuEfficiency::endJob() {
673  //  double efficiencySTA =(double)numberOfOverlappedStandAlone_/(double)numberOfMatchedZMuTrack_;
674  //  double errorEff_STA = sqrt( efficiencySTA*(1 - efficiencySTA)/(double)numberOfMatchedZMuTrack_);
675
676  double myTrackEff = 2.*numberOfMatchedSelectedZMuMu_/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuSta_);
677  double myErrTrackEff = sqrt(myTrackEff*(1-myTrackEff)/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuSta_));
678
679  double myStaEff = 2.*numberOfMatchedSelectedZMuMu_/(2.*numberOfMatchedSelectedZMuMu_+(double)numberOfMatchedSelectedZMuTrack_exclusive);
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  
720#include "FWCore/Framework/interface/MakerMacros.h"
721
722DEFINE_FWK_MODULE(ZMuMuEfficiency);
723