PageRenderTime 102ms CodeModel.GetById 19ms app.highlight 75ms RepoModel.GetById 1ms app.codeStats 1ms

/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMu_efficiencyAnalyzer.cc

https://github.com/aivanov-cern/cmssw
C++ | 748 lines | 603 code | 92 blank | 53 comment | 199 complexity | 6f0e80ffb8853dab62778ac98ba06ef8 MD5 | raw file
  1/* \class ZMuMu_efficieyAnalyzer
  2 * 
  3 * author: Davide Piccolo
  4 *
  5 * ZMuMu efficiency analyzer:
  6 * check muon reco efficiencies from MC truth, 
  7 *
  8 */
  9
 10#include "DataFormats/Common/interface/AssociationVector.h"
 11#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 12#include "DataFormats/Candidate/interface/CandMatchMap.h" 
 13#include "FWCore/Framework/interface/EDAnalyzer.h"
 14#include "DataFormats/Candidate/interface/Particle.h"
 15#include "DataFormats/Candidate/interface/Candidate.h"
 16#include "DataFormats/Candidate/interface/CandidateFwd.h"
 17#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
 18#include "DataFormats/MuonReco/interface/Muon.h"
 19#include "FWCore/ParameterSet/interface/ParameterSet.h"
 20#include "FWCore/Framework/interface/Event.h"
 21#include "FWCore/Framework/interface/EventSetup.h"
 22#include "FWCore/Utilities/interface/InputTag.h"
 23#include "DataFormats/Candidate/interface/OverlapChecker.h"
 24#include "DataFormats/Math/interface/deltaR.h"
 25#include "DataFormats/PatCandidates/interface/Muon.h" 
 26#include "DataFormats/PatCandidates/interface/GenericParticle.h"
 27#include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h"
 28#include "DataFormats/Common/interface/AssociationVector.h"
 29#include "DataFormats/PatCandidates/interface/PATObject.h"
 30
 31#include "TH1.h"
 32#include "TH2.h"
 33#include "TH3.h"
 34#include <vector>
 35using namespace edm;
 36using namespace std;
 37using namespace reco;
 38
 39class ZMuMu_efficiencyAnalyzer : public edm::EDAnalyzer {
 40public:
 41  ZMuMu_efficiencyAnalyzer(const edm::ParameterSet& pset);
 42private:
 43  virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
 44  bool check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
 45  float getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
 46  float getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
 47  float getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
 48  Particle::LorentzVector getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2); 
 49  virtual void endJob() override;
 50
 51  edm::InputTag zMuMu_, zMuMuMatchMap_; 
 52  edm::InputTag zMuStandAlone_, zMuStandAloneMatchMap_;
 53  edm::InputTag zMuTrack_, zMuTrackMatchMap_; 
 54  edm::InputTag muons_, muonMatchMap_, muonIso_;
 55  edm::InputTag tracks_, trackIso_;
 56  edm::InputTag genParticles_, primaryVertices_;
 57
 58  bool bothMuons_;
 59
 60  double etamax_, ptmin_, massMin_, massMax_, isoMax_;
 61
 62  // binning of entries array (at moment defined by hand and not in cfg file)
 63  unsigned int etaBins;
 64  unsigned int ptBins;
 65  double  etaRange[7];
 66  double  ptRange[5];
 67
 68  reco::CandidateBaseRef globalMuonCandRef_, trackMuonCandRef_, standAloneMuonCandRef_;
 69  OverlapChecker overlap_;
 70
 71  // general histograms
 72  TH1D *h_zmm_mass, *h_zmm2HLT_mass; 
 73  TH1D *h_zmm1HLTplus_mass, *h_zmmNotIsoplus_mass, *h_zmsplus_mass, *h_zmtplus_mass;
 74  TH1D *h_zmm1HLTminus_mass, *h_zmmNotIsominus_mass, *h_zmsminus_mass, *h_zmtminus_mass;
 75
 76  // global counters
 77  int nGlobalMuonsMatched_passed;    // total number of global muons MC matched and passing cuts (and triggered)
 78
 79  vector<TH1D *>  hmumu2HLTplus_eta, hmumu1HLTplus_eta, hmustaplus_eta, hmutrackplus_eta, hmumuNotIsoplus_eta;
 80  vector<TH1D *>  hmumu2HLTplus_pt, hmumu1HLTplus_pt, hmustaplus_pt, hmutrackplus_pt, hmumuNotIsoplus_pt;
 81  vector<TH1D *>  hmumu2HLTminus_eta, hmumu1HLTminus_eta, hmustaminus_eta, hmutrackminus_eta, hmumuNotIsominus_eta;
 82  vector<TH1D *>  hmumu2HLTminus_pt, hmumu1HLTminus_pt, hmustaminus_pt, hmutrackminus_pt, hmumuNotIsominus_pt;
 83};
 84
 85#include "FWCore/ServiceRegistry/interface/Service.h"
 86#include "CommonTools/UtilAlgos/interface/TFileService.h"
 87#include "DataFormats/Common/interface/Handle.h"
 88#include "DataFormats/Candidate/interface/Particle.h"
 89#include "DataFormats/Candidate/interface/Candidate.h"
 90#include "DataFormats/Common/interface/ValueMap.h"
 91#include "DataFormats/Candidate/interface/CandAssociation.h"
 92#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 93#include "DataFormats/Math/interface/LorentzVector.h"
 94#include "DataFormats/TrackReco/interface/Track.h"
 95#include <iostream>
 96#include <iterator>
 97#include <cmath>
 98using namespace std;
 99using namespace reco;
100using namespace edm;
101
102
103typedef edm::ValueMap<float> IsolationCollection;
104
105ZMuMu_efficiencyAnalyzer::ZMuMu_efficiencyAnalyzer(const ParameterSet& pset) : 
106  zMuMu_(pset.getParameter<InputTag>("zMuMu")), 
107  zMuMuMatchMap_(pset.getParameter<InputTag>("zMuMuMatchMap")), 
108  zMuStandAlone_(pset.getParameter<InputTag>("zMuStandAlone")), 
109  zMuStandAloneMatchMap_(pset.getParameter<InputTag>("zMuStandAloneMatchMap")), 
110  zMuTrack_(pset.getParameter<InputTag>("zMuTrack")), 
111  zMuTrackMatchMap_(pset.getParameter<InputTag>("zMuTrackMatchMap")), 
112  muons_(pset.getParameter<InputTag>("muons")), 
113  tracks_(pset.getParameter<InputTag>("tracks")), 
114  genParticles_(pset.getParameter<InputTag>( "genParticles" ) ),
115  primaryVertices_(pset.getParameter<InputTag>( "primaryVertices" ) ),
116
117  bothMuons_(pset.getParameter<bool>("bothMuons")), 
118
119  etamax_(pset.getUntrackedParameter<double>("etamax")),  
120  ptmin_(pset.getUntrackedParameter<double>("ptmin")), 
121  massMin_(pset.getUntrackedParameter<double>("zMassMin")), 
122  massMax_(pset.getUntrackedParameter<double>("zMassMax")), 
123  isoMax_(pset.getUntrackedParameter<double>("isomax")) { 
124  Service<TFileService> fs;
125
126  // general histograms
127  h_zmm_mass  = fs->make<TH1D>("zmm_mass","zmumu mass",100,0.,200.);
128  h_zmm2HLT_mass  = fs->make<TH1D>("zmm2HLT_mass","zmumu 2HLT mass",100,0.,200.);
129  h_zmm1HLTplus_mass  = fs->make<TH1D>("zmm1HLTplus_mass","zmumu 1HLT plus mass",100,0.,200.);
130  h_zmmNotIsoplus_mass  = fs->make<TH1D>("zmmNotIsoplus_mass","zmumu a least One Not Iso plus mass",100,0.,200.);
131  h_zmsplus_mass  = fs->make<TH1D>("zmsplus_mass","zmusta plus mass",100,0.,200.);
132  h_zmtplus_mass  = fs->make<TH1D>("zmtplus_mass","zmutrack plus mass",100,0.,200.);
133  h_zmm1HLTminus_mass  = fs->make<TH1D>("zmm1HLTminus_mass","zmumu 1HLT minus mass",100,0.,200.);
134  h_zmmNotIsominus_mass  = fs->make<TH1D>("zmmNotIsominus_mass","zmumu a least One Not Iso minus mass",100,0.,200.);
135  h_zmsminus_mass  = fs->make<TH1D>("zmsminus_mass","zmusta minus mass",100,0.,200.);
136  h_zmtminus_mass  = fs->make<TH1D>("zmtminus_mass","zmutrack minus mass",100,0.,200.);
137
138  cout << "primo" << endl;
139  // creating histograms for each Pt, eta interval
140
141  TFileDirectory etaDirectory = fs->mkdir("etaIntervals");   // in this directory will be saved all the histos of different eta intervals
142  TFileDirectory ptDirectory = fs->mkdir("ptIntervals");   // in this directory will be saved all the histos of different pt intervals
143
144  // binning of entries array (at moment defined by hand and not in cfg file)
145  etaBins = 6;
146  ptBins = 4;
147  double  etaRangeTmp[7] = {-2.,-1.2,-0.8,0.,0.8,1.2,2.};
148  double  ptRangeTmp[5] = {20.,40.,60.,80.,100.};
149  for (unsigned int i=0;i<=etaBins;i++) etaRange[i] = etaRangeTmp[i];
150  for (unsigned int i=0;i<=ptBins;i++) ptRange[i] = ptRangeTmp[i];
151
152  // eta histograms creation
153  cout << "eta istograms creation " << endl;
154
155  for (unsigned int i=0;i<etaBins;i++) {
156    cout << " bin eta plus  " << i << endl;
157    // muon plus
158    double range0 = etaRange[i];
159    double range1= etaRange[i+1];
160    char ap[30], bp[50];
161    sprintf(ap,"zmumu2HLTplus_etaRange%d",i);
162    sprintf(bp,"zmumu2HLT plus mass eta Range %f to %f",range0,range1);
163    cout << ap << "   " << bp << endl;
164    hmumu2HLTplus_eta.push_back(etaDirectory.make<TH1D>(ap,bp,200,0.,200.));
165    sprintf(ap,"zmumu1HLTplus_etaRange%d",i);
166    sprintf(bp,"zmumu1HLT plus mass eta Range %f to %f",range0,range1);
167    cout << ap << "   " << bp << endl;
168    hmumu1HLTplus_eta.push_back(etaDirectory.make<TH1D>(ap,bp,200,0.,200.));
169    sprintf(ap,"zmustaplus_etaRange%d",i);
170    sprintf(bp,"zmusta plus mass eta Range %f to %f",range0,range1);
171    cout << ap << "   " << bp << endl;
172    hmustaplus_eta.push_back(etaDirectory.make<TH1D>(ap,bp,50,0.,200.));
173    sprintf(ap,"zmutrackplus_etaRange%d",i);
174    sprintf(bp,"zmutrack plus mass eta Range %f to %f",range0,range1);
175    cout << ap << "   " << bp << endl;
176    hmutrackplus_eta.push_back(etaDirectory.make<TH1D>(ap,bp,100,0.,200.));
177    sprintf(ap,"zmumuNotIsoplus_etaRange%d",i);
178    sprintf(bp,"zmumuNotIso plus mass eta Range %f to %f",range0,range1);
179    cout << ap << "   " << bp << endl;
180    hmumuNotIsoplus_eta.push_back(etaDirectory.make<TH1D>(ap,bp,100,0.,200.));
181    // muon minus 
182    cout << " bin eta minus  " << i << endl;
183    char am[30], bm[50];
184    sprintf(am,"zmumu2HLTminus_etaRange%d",i);
185    sprintf(bm,"zmumu2HLT minus mass eta Range %f to %f",range0,range1);
186    cout << am << "   " << bm << endl;
187    hmumu2HLTminus_eta.push_back(etaDirectory.make<TH1D>(am,bm,200,0.,200.));
188    sprintf(am,"zmumu1HLTminus_etaRange%d",i);
189    sprintf(bm,"zmumu1HLT minus mass eta Range %f to %f",range0,range1);
190    cout << am << "   " << bm << endl;
191    hmumu1HLTminus_eta.push_back(etaDirectory.make<TH1D>(am,bm,200,0.,200.));
192    sprintf(am,"zmustaminus_etaRange%d",i);
193    sprintf(bm,"zmusta minus mass eta Range %f to %f",range0,range1);
194    cout << am << "   " << bm << endl;
195    hmustaminus_eta.push_back(etaDirectory.make<TH1D>(am,bm,50,0.,200.));
196    sprintf(am,"zmutrackminus_etaRange%d",i);
197    sprintf(bm,"zmutrack minus mass eta Range %f to %f",range0,range1);
198    cout << am << "   " << bm << endl;
199    hmutrackminus_eta.push_back(etaDirectory.make<TH1D>(am,bm,100,0.,200.));
200    sprintf(am,"zmumuNotIsominus_etaRange%d",i);
201    sprintf(bm,"zmumuNotIso minus mass eta Range %f to %f",range0,range1);
202    cout << am << "   " << bm << endl;
203    hmumuNotIsominus_eta.push_back(etaDirectory.make<TH1D>(am,bm,100,0.,200.));
204  } 
205
206  // pt histograms creation
207  cout << "pt istograms creation " << endl;
208
209  for (unsigned int i=0;i<ptBins;i++) {
210    double range0 = ptRange[i];
211    double range1= ptRange[i+1];
212    // muon plus
213    cout << " bin pt plus  " << i << endl;
214    char ap1[30], bp1[50];
215    sprintf(ap1,"zmumu2HLTplus_ptRange%d",i);
216    sprintf(bp1,"zmumu2HLT plus mass pt Range %f to %f",range0,range1);
217    cout << ap1 << "   " << bp1 << endl;
218    hmumu2HLTplus_pt.push_back(ptDirectory.make<TH1D>(ap1,bp1,200,0.,200.));
219    sprintf(ap1,"zmumu1HLTplus_ptRange%d",i);
220    sprintf(bp1,"zmumu1HLT plus mass pt Range %f to %f",range0,range1);
221    cout << ap1 << "   " << bp1 << endl;
222    hmumu1HLTplus_pt.push_back(ptDirectory.make<TH1D>(ap1,bp1,200,0.,200.));
223    sprintf(ap1,"zmustaplus_ptRange%d",i);
224    sprintf(bp1,"zmusta plus mass pt Range %f to %f",range0,range1);
225    cout << ap1 << "   " << bp1 << endl;
226    hmustaplus_pt.push_back(ptDirectory.make<TH1D>(ap1,bp1,50,0.,200.));
227    sprintf(ap1,"zmutrackplus_ptRange%d",i);
228    sprintf(bp1,"zmutrack plus mass pt Range %f to %f",range0,range1);
229    cout << ap1 << "   " << bp1 << endl;
230    hmutrackplus_pt.push_back(ptDirectory.make<TH1D>(ap1,bp1,100,0.,200.));
231    sprintf(ap1,"zmumuNotIsoplus_ptRange%d",i);
232    sprintf(bp1,"zmumuNotIso plus mass pt Range %f to %f",range0,range1);
233    cout << ap1 << "   " << bp1 << endl;
234    hmumuNotIsoplus_pt.push_back(ptDirectory.make<TH1D>(ap1,bp1,100,0.,200.));
235    // muon minus 
236    cout << " bin pt minus  " << i << endl;
237    char am1[30], bm1[50];
238    sprintf(am1,"zmumu2HLTminus_ptRange%d",i);
239    sprintf(bm1,"zmumu2HLT minus mass pt Range %f to %f",range0,range1);
240    cout << am1 << "   " << bm1 << endl;
241    hmumu2HLTminus_pt.push_back(ptDirectory.make<TH1D>(am1,bm1,200,0.,200.));
242    sprintf(am1,"zmumu1HLTminus_ptRange%d",i);
243    sprintf(bm1,"zmumu1HLT minus mass pt Range %f to %f",range0,range1);
244    cout << am1 << "   " << bm1 << endl;
245    hmumu1HLTminus_pt.push_back(ptDirectory.make<TH1D>(am1,bm1,200,0.,200.));
246    sprintf(am1,"zmustaminus_ptRange%d",i);
247    sprintf(bm1,"zmusta minus mass pt Range %f to %f",range0,range1);
248    cout << am1 << "   " << bm1 << endl;
249    hmustaminus_pt.push_back(ptDirectory.make<TH1D>(am1,bm1,50,0.,200.));
250    sprintf(am1,"zmutrackminus_ptRange%d",i);
251    sprintf(bm1,"zmutrack minus mass pt Range %f to %f",range0,range1);
252    cout << am1 << "   " << bm1 << endl;
253    hmutrackminus_pt.push_back(ptDirectory.make<TH1D>(am1,bm1,100,0.,200.));
254    sprintf(am1,"zmumuNotIsominus_ptRange%d",i);
255    sprintf(bm1,"zmumuNotIso minus mass pt Range %f to %f",range0,range1);
256    cout << am1 << "   " << bm1 << endl;
257    hmumuNotIsominus_pt.push_back(ptDirectory.make<TH1D>(am1,bm1,100,0.,200.));
258  } 
259
260  // clear global counters
261  nGlobalMuonsMatched_passed = 0;
262}
263
264void ZMuMu_efficiencyAnalyzer::analyze(const Event& event, const EventSetup& setup) {
265  Handle<CandidateView> zMuMu;  
266  Handle<GenParticleMatch> zMuMuMatchMap; //Map of Z made by Mu global + Mu global 
267  Handle<CandidateView> zMuStandAlone;  
268  Handle<GenParticleMatch> zMuStandAloneMatchMap; //Map of Z made by Mu + StandAlone
269  Handle<CandidateView> zMuTrack;  
270  Handle<GenParticleMatch> zMuTrackMatchMap; //Map of Z made by Mu + Track
271  Handle<CandidateView> muons; //Collection of Muons
272  Handle<CandidateView> tracks; //Collection of Tracks
273
274  Handle<GenParticleCollection> genParticles;  // Collection of Generatd Particles
275  Handle<reco::VertexCollection> primaryVertices;  // Collection of primary Vertices
276  
277  event.getByLabel(zMuMu_, zMuMu); 
278  event.getByLabel(zMuStandAlone_, zMuStandAlone); 
279  event.getByLabel(zMuTrack_, zMuTrack); 
280  event.getByLabel(genParticles_, genParticles);
281  event.getByLabel(primaryVertices_, primaryVertices);
282  event.getByLabel(muons_, muons); 
283  event.getByLabel(tracks_, tracks); 
284
285  /*
286  cout << "*********  zMuMu         size : " << zMuMu->size() << endl;
287  cout << "*********  zMuStandAlone size : " << zMuStandAlone->size() << endl;
288  cout << "*********  zMuTrack      size : " << zMuTrack->size() << endl;
289  cout << "*********  muons         size : " << muons->size() << endl; 	    
290  cout << "*********  tracks        size : " << tracks->size() << endl;
291  cout << "*********  vertices      size : " << primaryVertices->size() << endl;
292  */
293
294  //      std::cout<<"Run-> "<<event.id().run()<<std::endl;
295  //      std::cout<<"Event-> "<<event.id().event()<<std::endl; 
296
297
298
299  bool zMuMu_found = false;
300  // loop on ZMuMu
301  if (zMuMu->size() > 0 ) {
302    for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates
303      const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate
304      CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i);
305
306      const Candidate * lep0 = zMuMuCand.daughter( 0 );
307      const Candidate * lep1 = zMuMuCand.daughter( 1 );
308      const pat::Muon & muonDau0 = dynamic_cast<const pat::Muon &>(*lep0->masterClone());
309      double trkiso0 = muonDau0.trackIso();
310      const pat::Muon & muonDau1 = dynamic_cast<const pat::Muon &>(*lep1->masterClone());
311      double trkiso1 = muonDau1.trackIso();
312
313      // kinemtic variables
314      double pt0 = zMuMuCand.daughter(0)->pt();
315      double pt1 = zMuMuCand.daughter(1)->pt();
316      double eta0 = zMuMuCand.daughter(0)->eta();
317      double eta1 = zMuMuCand.daughter(1)->eta();
318      double charge0 = zMuMuCand.daughter(0)->charge();
319      double charge1 = zMuMuCand.daughter(1)->charge();
320      double mass = zMuMuCand.mass();
321
322      // HLT match
323      const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
324	muonDau0.triggerObjectMatchesByPath( "HLT_Mu9" );
325      const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 
326	muonDau1.triggerObjectMatchesByPath( "HLT_Mu9" );
327
328      bool trig0found = false;
329      bool trig1found = false;
330      if( mu0HLTMatches.size()>0 )
331	trig0found = true;
332      if( mu1HLTMatches.size()>0 )
333	trig1found = true;
334      
335      // kinematic selection
336
337      bool checkOppositeCharge = false;
338      if (charge0 != charge1) checkOppositeCharge = true;
339      if (pt0>ptmin_ && pt1>ptmin_ && abs(eta0)<etamax_ && abs(eta1)<etamax_ && mass>massMin_ && mass<massMax_ && checkOppositeCharge) {
340	if (trig0found || trig1found) { // at least one muon match HLT 
341	  zMuMu_found = true;           // Z found as global-global (so don't check Zms and Zmt)
342	  if (trkiso0 < isoMax_ && trkiso1 < isoMax_) { // both muons are isolated
343	    if (trig0found && trig1found) {  
344	      
345	      // ******************** category zmm 2 HLT ****************
346	      
347	      h_zmm2HLT_mass->Fill(mass);  
348	      h_zmm_mass->Fill(mass);  
349	      
350	      // check the cynematics to fill correct histograms
351	      
352	      for (unsigned int j=0;j<etaBins;j++) {  // eta Bins loop
353		double range0 = etaRange[j];
354		double range1= etaRange[j+1];
355
356		// eta histograms
357		
358		if (eta0>=range0 && eta0<range1)
359		  {
360		    if (charge0<0) hmumu2HLTminus_eta[j]->Fill(mass);  // mu- in bin eta
361		    if (charge0>0) hmumu2HLTplus_eta[j]->Fill(mass);  // mu+ in bin eta
362		  }
363		if (eta1>=range0 && eta1<range1)
364		  {
365		    if (charge1<0) hmumu2HLTminus_eta[j]->Fill(mass);  // mu- in bin eta
366		    if (charge1>0) hmumu2HLTplus_eta[j]->Fill(mass);  // mu+ in bin eta
367		  }
368	      } // end loop etaBins
369	      
370	      for (unsigned int j=0;j<ptBins;j++) {  // pt Bins loop
371		double range0pt = ptRange[j];
372		double range1pt = ptRange[j+1];
373		// pt histograms
374		if (pt0>=range0pt && pt0<range1pt)
375		  {
376		    if (charge0<0) hmumu2HLTminus_pt[j]->Fill(mass);  // mu- in bin eta
377		    if (charge0>0) hmumu2HLTplus_pt[j]->Fill(mass);  // mu+ in bin eta
378		  }
379		if (pt1>=range0pt && pt1<range1pt)
380		  {
381		    if (charge1<0) hmumu2HLTminus_pt[j]->Fill(mass);  // mu- in bin eta
382		    if (charge1>0) hmumu2HLTplus_pt[j]->Fill(mass);  // mu+ in bin eta
383		  }
384	      } // end loop  ptBins
385	      
386	    }  // ******************* end category zmm 2 HLT ****************
387	    
388	    if (!trig0found || !trig1found) { 
389	      // ****************** category zmm 1 HLT ****************** 
390	      h_zmm_mass->Fill(mass);  
391	      double eta = 9999;
392	      double pt = 9999;
393	      double charge = 0;
394	      if (trig0found) {
395		eta = eta1;       // check  muon not HLT matched
396		pt = pt1;
397		charge = charge1;
398	      } else {
399		eta = eta0;       
400		pt =pt0;
401		charge = charge0;
402	      }
403	      if (charge<0) h_zmm1HLTminus_mass->Fill(mass);  
404	      if (charge>0) h_zmm1HLTplus_mass->Fill(mass);  
405
406	      for (unsigned int j=0;j<etaBins;j++) {  // eta Bins loop
407		double range0 = etaRange[j];
408		double range1= etaRange[j+1];
409		// eta histograms fill the bin of the muon not HLT matched
410		if (eta>=range0 && eta<range1) 
411		  {
412		    if (charge<0) hmumu1HLTminus_eta[j]->Fill(mass);
413		    if (charge>0) hmumu1HLTplus_eta[j]->Fill(mass);
414		  }
415	      } // end loop etaBins
416	      for (unsigned int j=0;j<ptBins;j++) {  // pt Bins loop
417		double range0 = ptRange[j];
418		double range1= ptRange[j+1];
419		// pt histograms
420		if (pt>=range0 && pt<range1)
421		  {
422		    if (charge<0) hmumu1HLTminus_pt[j]->Fill(mass);
423		    if (charge>0) hmumu1HLTplus_pt[j]->Fill(mass);
424		  }
425	      } // end loop ptBins
426
427	    } // ****************** end category zmm 1 HLT ***************
428
429	  } else {  // one or both muons are not isolated
430	    // *****************  category zmumuNotIso **************** (per ora non studio iso vs eta e pt da capire meglio)
431
432	  } // end if both muons isolated
433	  
434	} // end if at least 1 HLT trigger found
435      }  // end if kinematic selection 
436
437
438    }  // end loop on ZMuMu cand
439  }    // end if ZMuMu size > 0
440
441  // loop on ZMuSta
442  bool zMuSta_found = false;
443  if (!zMuMu_found && zMuStandAlone->size() > 0 ) {
444    event.getByLabel(zMuStandAloneMatchMap_, zMuStandAloneMatchMap); 
445    for(unsigned int i = 0; i < zMuStandAlone->size(); ++i) { //loop on candidates
446      const Candidate & zMuStandAloneCand = (*zMuStandAlone)[i]; //the candidate
447      CandidateBaseRef zMuStandAloneCandRef = zMuStandAlone->refAt(i);
448      GenParticleRef zMuStandAloneMatch = (*zMuStandAloneMatchMap)[zMuStandAloneCandRef];
449
450      const Candidate * lep0 = zMuStandAloneCand.daughter( 0 );
451      const Candidate * lep1 = zMuStandAloneCand.daughter( 1 );
452      const pat::Muon & muonDau0 = dynamic_cast<const pat::Muon &>(*lep0->masterClone());
453      double trkiso0 = muonDau0.trackIso();
454      const pat::Muon & muonDau1 = dynamic_cast<const pat::Muon &>(*lep1->masterClone());
455      double trkiso1 = muonDau1.trackIso();
456      double pt0 = zMuStandAloneCand.daughter(0)->pt();
457      double pt1 = zMuStandAloneCand.daughter(1)->pt();
458      double eta0 = zMuStandAloneCand.daughter(0)->eta();
459      double eta1 = zMuStandAloneCand.daughter(1)->eta();
460      double charge0 = zMuStandAloneCand.daughter(0)->charge();
461      double charge1 = zMuStandAloneCand.daughter(1)->charge();
462      double mass = zMuStandAloneCand.mass();
463
464      // HLT match
465      const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
466	muonDau0.triggerObjectMatchesByPath( "HLT_Mu9" );
467      const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 
468	muonDau1.triggerObjectMatchesByPath( "HLT_Mu9" );
469
470      bool trig0found = false;
471      bool trig1found = false;
472      if( mu0HLTMatches.size()>0 )
473	trig0found = true;
474      if( mu1HLTMatches.size()>0 )
475	trig1found = true;
476
477      // check HLT match of Global muon and save eta, pt of second muon (standAlone)
478      bool trigGlbfound = false;
479      double pt =999.;
480      double eta = 999.;
481      double charge = 0;
482      if (muonDau0.isGlobalMuon()) {
483	trigGlbfound = trig0found;
484	pt = pt1;
485	eta = eta1;
486	charge = charge1;
487      }
488      if (muonDau1.isGlobalMuon()) {
489	trigGlbfound = trig1found;
490	pt = pt0;
491	eta = eta0;
492	charge = charge0;
493      }
494
495      bool checkOppositeCharge = false;
496      if (charge0 != charge1) checkOppositeCharge = true;
497
498      if (checkOppositeCharge && trigGlbfound && pt0>ptmin_ && pt1>ptmin_ && abs(eta0)<etamax_ && abs(eta1)<etamax_ && mass>massMin_ && mass<massMax_ && trkiso0<isoMax_ && trkiso1<isoMax_ ) {  // global mu match HLT + kinematic cuts + opposite charge
499
500	if (charge<0) h_zmsminus_mass->Fill(mass);
501	if (charge>0) h_zmsplus_mass->Fill(mass);
502
503	for (unsigned int j=0;j<etaBins;j++) {  // eta Bins loop
504	  double range0 = etaRange[j];
505	  double range1= etaRange[j+1];
506	  // eta histograms
507	  if (eta>=range0 && eta<range1) {
508	    if (charge<0)  hmustaminus_eta[j]->Fill(mass);
509	    if (charge>0)  hmustaplus_eta[j]->Fill(mass);
510	  }
511	} // end loop etaBins
512	for (unsigned int j=0;j<ptBins;j++) {  // pt Bins loop
513	  double range0 = ptRange[j];
514	  double range1= ptRange[j+1];
515	  // pt histograms
516	  if (pt>=range0 && pt<range1) {
517	    if (charge<0)  hmustaminus_pt[j]->Fill(mass);
518	    if (charge>0)  hmustaplus_pt[j]->Fill(mass);
519	  }
520	} // end loop ptBins
521	
522      } // end if trigGlbfound + kinecuts + OppostieCharge
523    }  // end loop on ZMuStandAlone cand
524  }    // end if ZMuStandAlone size > 0
525
526
527  // loop on ZMuTrack
528  //  bool zMuTrack_found = false;
529  if (!zMuMu_found && !zMuSta_found && zMuTrack->size() > 0 ) {
530    event.getByLabel(zMuTrackMatchMap_, zMuTrackMatchMap); 
531    for(unsigned int i = 0; i < zMuTrack->size(); ++i) { //loop on candidates
532      const Candidate & zMuTrackCand = (*zMuTrack)[i]; //the candidate
533      CandidateBaseRef zMuTrackCandRef = zMuTrack->refAt(i);
534      const Candidate * lep0 = zMuTrackCand.daughter( 0 );
535      const Candidate * lep1 = zMuTrackCand.daughter( 1 );
536      const pat::Muon & muonDau0 = dynamic_cast<const pat::Muon &>(*lep0->masterClone());
537      double trkiso0 = muonDau0.trackIso();
538      const pat::GenericParticle & trackDau1 = dynamic_cast<const pat::GenericParticle &>(*lep1->masterClone());
539      double trkiso1 = trackDau1.trackIso();
540      double pt0 = zMuTrackCand.daughter(0)->pt();
541      double pt1 = zMuTrackCand.daughter(1)->pt();
542      double eta0 = zMuTrackCand.daughter(0)->eta();
543      double eta1 = zMuTrackCand.daughter(1)->eta();
544      double charge0 = zMuTrackCand.daughter(0)->charge();
545      double charge1 = zMuTrackCand.daughter(1)->charge();
546      double mass = zMuTrackCand.mass();
547
548      // HLT match (check just dau0 the global)
549      const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
550	muonDau0.triggerObjectMatchesByPath( "HLT_Mu9" );
551
552      bool trig0found = false;
553      if( mu0HLTMatches.size()>0 )
554	trig0found = true;
555
556      bool checkOppositeCharge = false;
557      if (charge0 != charge1) checkOppositeCharge = true;
558
559      if (checkOppositeCharge && trig0found && pt0>ptmin_ && pt1>ptmin_ && abs(eta0)<etamax_ && abs(eta1)<etamax_ && mass>massMin_ && mass<massMax_ && trkiso0<isoMax_ && trkiso1<isoMax_ ) {  // global mu match HLT + kinematic cuts + opposite charge
560
561	if (charge1<0) h_zmtminus_mass->Fill(mass);
562	if (charge1>0) h_zmtplus_mass->Fill(mass);
563
564	for (unsigned int j=0;j<etaBins;j++) {  // eta Bins loop
565	  double range0 = etaRange[j];
566	  double range1= etaRange[j+1];
567	  // eta histograms
568	  if (eta1>=range0 && eta1<range1) {
569	    if (charge1<0)  hmutrackminus_eta[j]->Fill(mass);  // just check muon1 (mu0 is global by definition)
570	    if (charge1>0)  hmutrackplus_eta[j]->Fill(mass);  // just check muon1 (mu0 is global by definition)
571	  }
572	} // end loop etaBins
573	for (unsigned int j=0;j<ptBins;j++) {  // pt Bins loop
574	  double range0 = ptRange[j];
575	  double range1= ptRange[j+1];
576	  // pt histograms
577	  if (pt1>=range0 && pt1<range1) {
578	    if (charge1<0)  hmutrackminus_pt[j]->Fill(mass);  // just check muon1 (mu0 is global by definition)
579	    if (charge1>0)  hmutrackplus_pt[j]->Fill(mass);  // just check muon1 (mu0 is global by definition)
580	  }
581	} // end loop ptBins
582	
583      } // end if trig0found
584
585      
586    }  // end loop on ZMuTrack cand
587  }    // end if ZMuTrack size > 0
588
589}       // end analyze
590
591bool ZMuMu_efficiencyAnalyzer::check_ifZmumu(const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
592{
593  int partId0 = dauGen0->pdgId();
594  int partId1 = dauGen1->pdgId();
595  int partId2 = dauGen2->pdgId();
596  bool muplusFound=false;
597  bool muminusFound=false;
598  bool ZFound=false;
599  if (partId0==13 || partId1==13 || partId2==13) muminusFound=true;
600  if (partId0==-13 || partId1==-13 || partId2==-13) muplusFound=true;
601  if (partId0==23 || partId1==23 || partId2==23) ZFound=true;
602  return muplusFound*muminusFound*ZFound;   
603}
604 
605float ZMuMu_efficiencyAnalyzer::getParticlePt(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
606{
607  int partId0 = dauGen0->pdgId();
608  int partId1 = dauGen1->pdgId();
609  int partId2 = dauGen2->pdgId();
610  float ptpart=0.;
611  if (partId0 == ipart) {
612    for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
613      const Candidate * dauMuGen = dauGen0->daughter(k);
614      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
615	ptpart = dauMuGen->pt();
616      }
617    }
618  }
619  if (partId1 == ipart) {
620    for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
621      const Candidate * dauMuGen = dauGen1->daughter(k);
622      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
623	ptpart = dauMuGen->pt();
624      }
625    }
626  }
627  if (partId2 == ipart) {
628    for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
629      const Candidate * dauMuGen = dauGen2->daughter(k);
630      if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
631	ptpart = dauMuGen->pt();
632      }
633    }
634  }
635  return ptpart;
636}
637 
638float ZMuMu_efficiencyAnalyzer::getParticleEta(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
639{
640  int partId0 = dauGen0->pdgId();
641  int partId1 = dauGen1->pdgId();
642  int partId2 = dauGen2->pdgId();
643  float etapart=0.;
644  if (partId0 == ipart) {
645    for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
646      const Candidate * dauMuGen = dauGen0->daughter(k);
647      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
648	etapart = dauMuGen->eta();
649      }
650    }
651  }
652  if (partId1 == ipart) {
653    for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
654      const Candidate * dauMuGen = dauGen1->daughter(k);
655      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
656	etapart = dauMuGen->eta();
657      }
658    }
659  }
660  if (partId2 == ipart) {
661    for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
662      const Candidate * dauMuGen = dauGen2->daughter(k);
663      if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
664	etapart = dauMuGen->eta();
665      }
666    }
667  }
668  return etapart;
669}
670
671float ZMuMu_efficiencyAnalyzer::getParticlePhi(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
672{
673  int partId0 = dauGen0->pdgId();
674  int partId1 = dauGen1->pdgId();
675  int partId2 = dauGen2->pdgId();
676  float phipart=0.;
677  if (partId0 == ipart) {
678    for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
679      const Candidate * dauMuGen = dauGen0->daughter(k);
680      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
681	phipart = dauMuGen->phi();
682      }
683    }
684  }
685  if (partId1 == ipart) {
686    for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
687      const Candidate * dauMuGen = dauGen1->daughter(k);
688      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
689	phipart = dauMuGen->phi();
690      }
691    }
692  }
693  if (partId2 == ipart) {
694    for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
695      const Candidate * dauMuGen = dauGen2->daughter(k);
696      if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
697	phipart = dauMuGen->phi();
698      }
699    }
700  }
701  return phipart;
702}
703
704Particle::LorentzVector ZMuMu_efficiencyAnalyzer::getParticleP4(const int ipart, const Candidate * dauGen0, const Candidate * dauGen1, const Candidate * dauGen2)
705{
706  int partId0 = dauGen0->pdgId();
707  int partId1 = dauGen1->pdgId();
708  int partId2 = dauGen2->pdgId();
709  Particle::LorentzVector p4part(0.,0.,0.,0.);
710  if (partId0 == ipart) {
711    for(unsigned int k = 0; k < dauGen0->numberOfDaughters(); ++k) {
712      const Candidate * dauMuGen = dauGen0->daughter(k);
713      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
714	p4part = dauMuGen->p4();
715      }
716    }
717  }
718  if (partId1 == ipart) {
719    for(unsigned int k = 0; k < dauGen1->numberOfDaughters(); ++k) {
720      const Candidate * dauMuGen = dauGen1->daughter(k);
721      if(dauMuGen->pdgId() == ipart && dauMuGen->status() ==1) {
722	p4part = dauMuGen->p4();
723      }
724    }
725  }
726  if (partId2 == ipart) {
727    for(unsigned int k = 0; k < dauGen2->numberOfDaughters(); ++k) {
728      const Candidate * dauMuGen = dauGen2->daughter(k);
729      if(abs(dauMuGen->pdgId()) == ipart && dauMuGen->status() ==1) {
730	p4part = dauMuGen->p4();
731      }
732    }
733  }
734  return p4part;
735}
736 
737
738
739void ZMuMu_efficiencyAnalyzer::endJob() {
740  
741  
742 
743}
744  
745#include "FWCore/Framework/interface/MakerMacros.h"
746
747DEFINE_FWK_MODULE(ZMuMu_efficiencyAnalyzer);
748