PageRenderTime 82ms CodeModel.GetById 12ms app.highlight 64ms RepoModel.GetById 1ms app.codeStats 0ms

/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMu_radiative_analysis.cc

https://github.com/aivanov-cern/cmssw
C++ | 558 lines | 476 code | 59 blank | 23 comment | 79 complexity | cd7c36f24c1efd326ae42582bcc49626 MD5 | raw file
  1/* \class ZMuMu_Radiative_analyzer
  2 * 
  3 * author: Pasquale Noli
  4 *
  5 * ZMuMu Radiative analyzer
  6 * 
  7 *
  8 */
  9#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 10#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
 11#include "DataFormats/MuonReco/interface/Muon.h"
 12#include "DataFormats/Common/interface/AssociationVector.h"
 13#include "DataFormats/Candidate/interface/CandMatchMap.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/Candidate/interface/OverlapChecker.h"
 18#include "DataFormats/PatCandidates/interface/Muon.h" 
 19#include "DataFormats/PatCandidates/interface/Lepton.h" 
 20#include "DataFormats/PatCandidates/interface/GenericParticle.h"
 21#include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h"
 22#include "DataFormats/PatCandidates/interface/PATObject.h"
 23#include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
 24#include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
 25#include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
 26#include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
 27#include "DataFormats/PatCandidates/interface/Isolation.h"
 28#include "DataFormats/Common/interface/Handle.h"
 29#include "DataFormats/Common/interface/ValueMap.h"
 30#include "DataFormats/Candidate/interface/CandAssociation.h"
 31#include "FWCore/Framework/interface/EDAnalyzer.h"
 32#include "FWCore/ParameterSet/interface/ParameterSet.h"
 33#include "FWCore/Framework/interface/Event.h"
 34#include "FWCore/Framework/interface/EventSetup.h"
 35#include "FWCore/Utilities/interface/InputTag.h"
 36#include "FWCore/ServiceRegistry/interface/Service.h"
 37#include "CommonTools/UtilAlgos/interface/TFileService.h"
 38#include <iostream>
 39#include <iterator>
 40#include <cmath>
 41#include <vector>
 42#include "TH1.h"
 43#include "TH2.h"
 44#include "TH3.h"
 45
 46
 47using namespace edm;
 48using namespace std;
 49using namespace reco;
 50using namespace isodeposit;
 51
 52class ZMuMu_Radiative_analyzer : public edm::EDAnalyzer {
 53public:
 54  ZMuMu_Radiative_analyzer(const edm::ParameterSet& pset);
 55private:
 56  virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
 57  virtual void endJob() override;
 58
 59  edm::InputTag zMuMu_, zMuMuMatchMap_, zMuTk_, zMuTkMatchMap_,zMuSa_, zMuSaMatchMap_;
 60  double dRVeto_, dRTrk_, ptThreshold_;
 61  //histograms 
 62  TH1D *h_zmass_FSR,*h_zmass_no_FSR;
 63  TH1D *h_zMuSamass_FSR,*h_zMuSamass_no_FSR;
 64  TH1D *h_zMuTkmass_FSR,*h_zMuTkmass_no_FSR;
 65  TH1D *h_Iso_,*h_Iso_FSR_ ;
 66  TH3D *h_Iso_3D_, *h_Iso_FSR_3D_; 
 67  TH2D *h_staProbe_pt_eta_no_FSR_, *h_staProbe_pt_eta_FSR_, *h_ProbeOk_pt_eta_no_FSR_, *h_ProbeOk_pt_eta_FSR_;
 68  TH1D *h_trackProbe_eta_no_FSR, *h_trackProbe_pt_no_FSR, *h_staProbe_eta_no_FSR, *h_staProbe_pt_no_FSR, *h_ProbeOk_eta_no_FSR, *h_ProbeOk_pt_no_FSR;
 69  TH1D *h_trackProbe_eta_FSR, *h_trackProbe_pt_FSR, *h_staProbe_eta_FSR, *h_staProbe_pt_FSR, *h_ProbeOk_eta_FSR, *h_ProbeOk_pt_FSR;
 70  //boolean 
 71  bool FSR_mu, FSR_tk, FSR_mu0, FSR_mu1;
 72  bool trig0found, trig1found;
 73  //counter
 74  int  zmmcounter , zmscounter, zmtcounter, evntcounter;
 75 };
 76
 77
 78typedef edm::ValueMap<float> IsolationCollection;
 79
 80ZMuMu_Radiative_analyzer::ZMuMu_Radiative_analyzer(const ParameterSet& pset) : 
 81  zMuMu_(pset.getParameter<InputTag>("zMuMu")), 
 82  zMuMuMatchMap_(pset.getParameter<InputTag>("zMuMuMatchMap")),
 83  zMuTk_(pset.getParameter<InputTag>("zMuTk")), 
 84  zMuTkMatchMap_(pset.getParameter<InputTag>("zMuTkMatchMap")),
 85  zMuSa_(pset.getParameter<InputTag>("zMuSa")), 
 86  zMuSaMatchMap_(pset.getParameter<InputTag>("zMuSaMatchMap")),
 87  dRVeto_(pset.getUntrackedParameter<double>("veto")),
 88  dRTrk_(pset.getUntrackedParameter<double>("deltaRTrk")),
 89  ptThreshold_(pset.getUntrackedParameter<double>("ptThreshold")){ 
 90  zmmcounter=0;
 91  zmscounter=0;
 92  zmtcounter=0;
 93  evntcounter=0;
 94  Service<TFileService> fs;
 95   
 96  // general histograms
 97  h_zmass_FSR= fs->make<TH1D>("h_zmass_FRS","Invariant Z mass distribution",200,0,200);
 98  h_zmass_no_FSR= fs->make<TH1D>("h_zmass_no_FSR","Invariant Z mass distribution",200,0,200);
 99  h_zMuSamass_FSR= fs->make<TH1D>("h_zMuSamass_FRS","Invariant Z mass distribution",200,0,200);
100  h_zMuSamass_no_FSR= fs->make<TH1D>("h_zMuSamass_no_FSR","Invariant Z mass distribution",200,0,200);
101  h_zMuTkmass_FSR= fs->make<TH1D>("h_zMuTkmass_FRS","Invariant Z mass distribution",200,0,200);
102  h_zMuTkmass_no_FSR= fs->make<TH1D>("h_zMuTkmass_no_FSR","Invariant Z mass distribution",200,0,200);
103  h_Iso_= fs->make<TH1D>("h_iso","Isolation distribution of muons without FSR",100,0,20);
104  h_Iso_FSR_= fs->make<TH1D>("h_iso_FSR","Isolation distribution of muons with FSR ",100,0,20);
105  h_Iso_3D_= fs->make<TH3D>("h_iso_3D","Isolation distribution of muons without FSR",100,20,100,100,-2.0,2.0,100,0,20);
106  h_Iso_FSR_3D_= fs->make<TH3D>("h_iso_FSR_3D","Isolation distribution of muons with FSR ",100,20,100,100,-2.0,2.0,100,0,20);
107  h_staProbe_pt_eta_no_FSR_= fs->make<TH2D>("h_staProbe_pt_eta_no_FSR","Pt vs Eta StandAlone without FSR ",100,20,100,100,-2.0,2.0);
108  h_staProbe_pt_eta_FSR_= fs->make<TH2D>("h_staProbe_pt_eta_FSR","Pt vs Eta StandAlone with FSR ",100,20,100,100,-2.0,2.0);
109  h_ProbeOk_pt_eta_no_FSR_= fs->make<TH2D>("h_ProbeOk_pt_eta_no_FSR","Pt vs Eta probeOk without FSR ",100,20,100,100,-2.0,2.0);
110  h_ProbeOk_pt_eta_FSR_= fs->make<TH2D>("h_ProbeOk_pt_eta_FSR","Pt vs Eta probeOk with FSR ",100,20,100,100,-2.0,2.0);
111  
112  h_trackProbe_eta_no_FSR = fs->make<TH1D>("trackProbeEta_no_FSR","Eta of tracks",100,-2.0,2.0);
113  h_trackProbe_pt_no_FSR = fs->make<TH1D>("trackProbePt_no_FSR","Pt of tracks",100,0,100);
114  h_staProbe_eta_no_FSR = fs->make<TH1D>("standAloneProbeEta_no_FSR","Eta of standAlone",100,-2.0,2.0);
115  h_staProbe_pt_no_FSR = fs->make<TH1D>("standAloneProbePt_no_FSR","Pt of standAlone",100,0,100);
116  h_ProbeOk_eta_no_FSR = fs->make<TH1D>("probeOkEta_no_FSR","Eta of probe Ok",100,-2.0,2.0);
117  h_ProbeOk_pt_no_FSR = fs->make<TH1D>("probeOkPt_no_FSR","Pt of probe ok",100,0,100);
118
119  h_trackProbe_eta_FSR = fs->make<TH1D>("trackProbeEta_FSR","Eta of tracks",100,-2.0,2.0);
120  h_trackProbe_pt_FSR  = fs->make<TH1D>("trackProbePt_FSR","Pt of tracks",100,0,100);
121  h_staProbe_eta_FSR  = fs->make<TH1D>("standAloneProbeEta_FSR","Eta of standAlone",100,-2.0,2.0);
122  h_staProbe_pt_FSR  = fs->make<TH1D>("standAloneProbePt_FSR","Pt of standAlone",100,0,100);
123  h_ProbeOk_eta_FSR  = fs->make<TH1D>("probeOkEta_FSR","Eta of probe Ok",100,-2.0,2.0);
124  h_ProbeOk_pt_FSR  = fs->make<TH1D>("probeOkPt_FSR","Pt of probe ok",100,0,100);
125
126
127
128}
129
130void ZMuMu_Radiative_analyzer::analyze(const Event& event, const EventSetup& setup) {
131  evntcounter++;
132  Handle<CandidateView> zMuMu;                //Collection of Z made by  Mu global + Mu global 
133  Handle<GenParticleMatch> zMuMuMatchMap;     //Map of Z made by Mu global + Mu global with MC
134  event.getByLabel(zMuMu_, zMuMu); 
135  Handle<CandidateView> zMuTk;                //Collection of Z made by  Mu global + Track 
136  Handle<GenParticleMatch> zMuTkMatchMap;
137  event.getByLabel(zMuTk_, zMuTk); 
138  Handle<CandidateView> zMuSa;                //Collection of Z made by  Mu global + Sa 
139  Handle<GenParticleMatch> zMuSaMatchMap;
140  event.getByLabel(zMuSa_, zMuSa); 
141  cout << "**********  New Event  ***********"<<endl; 
142  // ZMuMu
143  if (zMuMu->size() > 0 ) {
144    event.getByLabel(zMuMuMatchMap_, zMuMuMatchMap); 
145     for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates
146     
147      const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate
148      CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i);
149
150       
151      CandidateBaseRef dau0 = zMuMuCand.daughter(0)->masterClone();
152      CandidateBaseRef dau1 = zMuMuCand.daughter(1)->masterClone();
153      const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon
154      const pat::Muon& mu1 = dynamic_cast<const pat::Muon&>(*dau1);
155            
156      double zmass= zMuMuCand.mass();
157      double pt0 = mu0.pt();
158      double pt1 = mu1.pt();
159      double eta0 = mu0.eta();
160      double eta1 = mu1.eta();
161      if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){
162	GenParticleRef zMuMuMatch = (*zMuMuMatchMap)[zMuMuCandRef];
163	if(zMuMuMatch.isNonnull()) {  // ZMuMu matched
164	  zmmcounter++;
165	  cout<<"         Zmumu cuts && matched" <<endl;
166	  FSR_mu0 = false;
167	  FSR_mu1 = false;
168	  
169	  //Isodeposit
170	  const pat::IsoDeposit * mu0TrackIso =mu0.isoDeposit(pat::TrackIso);
171	  const pat::IsoDeposit * mu1TrackIso =mu1.isoDeposit(pat::TrackIso);
172	  Direction mu0Dir = Direction(mu0.eta(),mu0.phi());
173	  Direction mu1Dir = Direction(mu1.eta(),mu1.phi());
174	  
175	  reco::IsoDeposit::AbsVetos vetos_mu0;
176	  vetos_mu0.push_back(new ConeVeto( mu0Dir, dRVeto_ ));
177	  vetos_mu0.push_back(new ThresholdVeto( ptThreshold_ ));
178	  
179	  reco::IsoDeposit::AbsVetos vetos_mu1;
180	  vetos_mu1.push_back(new ConeVeto( mu1Dir, dRVeto_ ));
181	  vetos_mu1.push_back(new ThresholdVeto( ptThreshold_ ));
182	  
183	  double  Tracker_isovalue_mu0 = mu0TrackIso->sumWithin(dRTrk_,vetos_mu0);
184	  double  Tracker_isovalue_mu1 = mu1TrackIso->sumWithin(dRTrk_,vetos_mu1);
185	 
186	  //trigger study 
187	  trig0found = false;
188	  trig1found = false;
189	  const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
190	    mu0.triggerObjectMatchesByPath( "HLT_Mu9" );
191	  const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 
192	    mu1.triggerObjectMatchesByPath( "HLT_Mu9" );
193	  if( mu0HLTMatches.size()>0 )
194	    trig0found = true;
195	  if( mu1HLTMatches.size()>0 )
196	    trig1found = true;
197      
198	  //MonteCarlo Study
199	  const reco::GenParticle * muMc0 = mu0.genLepton();
200	  const reco::GenParticle * muMc1 = mu1.genLepton();
201	  const Candidate * motherMu0 =  muMc0->mother();
202	  int num_dau_muon0 = motherMu0->numberOfDaughters();
203	  const Candidate * motherMu1 =  muMc1->mother();
204	  int num_dau_muon1 = motherMu1->numberOfDaughters();
205	  cout<<"         muone0"<<endl;
206	  cout<<"         num di daughters = "<< num_dau_muon0 <<endl;
207	  if( num_dau_muon0 > 1 ){
208	    for(int j = 0; j <  num_dau_muon0; ++j){
209	      int id =motherMu0 ->daughter(j)->pdgId();
210	      cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
211	      if(id == 22) FSR_mu0=true;
212	    }
213	  }//end check of gamma
214
215	  cout<<"         muone1"<<endl;
216	  cout<<"         num di daughters = "<< num_dau_muon1 <<endl;
217	  if( num_dau_muon1 > 1 ){
218	    for(int j = 0; j <  num_dau_muon1; ++j){
219	      int id = motherMu1->daughter(j)->pdgId(); 
220	      cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
221	      if(id == 22) FSR_mu1=true;
222	    }
223	  }//end check of gamma
224	
225	  if(FSR_mu0 || FSR_mu1 )h_zmass_FSR->Fill(zmass);
226	  else h_zmass_no_FSR->Fill(zmass);
227
228	  if (trig1found) {       // check efficiency of muon0 not imposing the trigger on it 
229	    cout<<"muon 1 is triggered "<<endl;
230	    if(FSR_mu0){
231	      cout<< "and muon 0 does FSR"<<endl;
232	      h_trackProbe_eta_FSR->Fill(eta0);
233	      h_trackProbe_pt_FSR->Fill(pt0);
234	      h_staProbe_eta_FSR->Fill(eta0);
235	      h_staProbe_pt_FSR->Fill(pt0);
236	      h_staProbe_pt_eta_FSR_->Fill(pt0,eta0);
237	      h_ProbeOk_eta_FSR->Fill(eta0);
238	      h_ProbeOk_pt_FSR->Fill(pt0);
239	      h_ProbeOk_pt_eta_FSR_->Fill(pt0,eta0);
240	    }else{
241	      cout<<"and muon 0 doesn't FSR"<<endl;
242	      h_trackProbe_eta_no_FSR->Fill(eta0);
243	      h_trackProbe_pt_no_FSR->Fill(pt0);
244	      h_staProbe_eta_no_FSR->Fill(eta0);
245	      h_staProbe_pt_no_FSR->Fill(pt0);
246	      h_staProbe_pt_eta_no_FSR_->Fill(pt0,eta0);
247	      h_ProbeOk_eta_no_FSR->Fill(eta0);
248	      h_ProbeOk_pt_no_FSR->Fill(pt0);
249	      h_ProbeOk_pt_eta_no_FSR_->Fill(pt0,eta0);
250	    }
251	    if(FSR_mu0){
252	      h_Iso_FSR_->Fill(Tracker_isovalue_mu0);
253	      h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
254	    }
255	    else{
256	      h_Iso_->Fill(Tracker_isovalue_mu0);
257	      h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
258	    }
259	  }
260	  if (trig0found) {         // check efficiency of muon1 not imposing the trigger on it 
261	    cout<<"muon 0 is triggered"<<endl;
262	    if(FSR_mu1){ 
263	      cout<<"and muon 1  does FSR"<<endl;
264	      h_trackProbe_eta_FSR->Fill(eta1);
265	      h_staProbe_eta_FSR->Fill(eta1);
266	      h_trackProbe_pt_FSR->Fill(pt1);
267	      h_staProbe_pt_FSR->Fill(pt1);
268	      h_ProbeOk_eta_FSR->Fill(eta1);
269	      h_ProbeOk_pt_FSR->Fill(pt1);
270	      h_staProbe_pt_eta_FSR_->Fill(pt1,eta1);
271	      h_ProbeOk_pt_eta_FSR_->Fill(pt1,eta1);
272
273	    }else{
274	      cout<<"and muon 1 doesn't FSR"<<endl;
275	      h_trackProbe_eta_no_FSR->Fill(eta1);
276	      h_staProbe_eta_no_FSR->Fill(eta1);
277	      h_trackProbe_pt_no_FSR->Fill(pt1);
278	      h_staProbe_pt_no_FSR->Fill(pt1);
279	      h_ProbeOk_eta_no_FSR->Fill(eta1);
280	      h_ProbeOk_pt_no_FSR->Fill(pt1);
281	      h_staProbe_pt_eta_no_FSR_->Fill(pt1,eta1);
282	      h_ProbeOk_pt_eta_no_FSR_->Fill(pt1,eta1);
283
284
285	    }
286	    if(FSR_mu1){
287	      h_Iso_FSR_->Fill(Tracker_isovalue_mu1);
288	      h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
289	    }else{ 
290	      h_Iso_->Fill(Tracker_isovalue_mu1);
291	      h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
292	    }
293	  }
294	}// end MC match
295      }//end of cuts
296     }// end loop on ZMuMu cand
297  }// end if ZMuMu size > 0
298  
299  // ZMuSa
300  if (zMuSa->size() > 0 ) {
301    event.getByLabel(zMuSaMatchMap_, zMuSaMatchMap); 
302     for(unsigned int i = 0; i < zMuSa->size(); ++i) { //loop on candidates
303     
304      const Candidate & zMuSaCand = (*zMuSa)[i]; //the candidate
305      CandidateBaseRef zMuSaCandRef = zMuSa->refAt(i);
306      const Candidate *  lep0 =zMuSaCand.daughter(0);  
307      const Candidate *  lep1 =zMuSaCand.daughter(1);  
308      CandidateBaseRef dau0 = lep0->masterClone();
309      CandidateBaseRef dau1 = lep1->masterClone();
310      const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon
311      const pat::Muon& mu1 = dynamic_cast<const pat::Muon&>(*dau1);
312            
313      double zmass= zMuSaCand.mass();
314      double pt0 = mu0.pt();
315      double pt1 = mu1.pt();
316      double eta0 = mu0.eta();
317      double eta1 = mu1.eta();
318      if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){
319	GenParticleRef zMuSaMatch = (*zMuSaMatchMap)[zMuSaCandRef];
320	if(zMuSaMatch.isNonnull()) {  // ZMuSa matched
321	  cout<<"         Zmusa cuts && matched" <<endl;
322	  FSR_mu0 = false;
323	  FSR_mu1 = false;
324	  zmscounter++;	  
325	  //Isodeposit
326	  const pat::IsoDeposit * mu0TrackIso =mu0.isoDeposit(pat::TrackIso);
327	  const pat::IsoDeposit * mu1TrackIso =mu1.isoDeposit(pat::TrackIso);
328	  Direction mu0Dir = Direction(mu0.eta(),mu0.phi());
329	  Direction mu1Dir = Direction(mu1.eta(),mu1.phi());
330	  
331	  reco::IsoDeposit::AbsVetos vetos_mu0;
332	  vetos_mu0.push_back(new ConeVeto( mu0Dir, dRVeto_ ));
333	  vetos_mu0.push_back(new ThresholdVeto( ptThreshold_ ));
334	  
335	  reco::IsoDeposit::AbsVetos vetos_mu1;
336	  vetos_mu1.push_back(new ConeVeto( mu1Dir, dRVeto_ ));
337	  vetos_mu1.push_back(new ThresholdVeto( ptThreshold_ ));
338	  
339	  double  Tracker_isovalue_mu0 = mu0TrackIso->sumWithin(dRTrk_,vetos_mu0);
340	  double  Tracker_isovalue_mu1 = mu1TrackIso->sumWithin(dRTrk_,vetos_mu1);
341
342	  // HLT match (check just dau0 the global)
343	  const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
344	    mu0.triggerObjectMatchesByPath( "HLT_Mu9" );
345	  const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 
346	    mu1.triggerObjectMatchesByPath( "HLT_Mu9" );
347	  trig0found = false;
348	  trig1found = false;
349	  if( mu0HLTMatches.size()>0 )
350	    trig0found = true;
351	  if( mu1HLTMatches.size()>0 )
352	    trig1found = true;
353	  
354	  //MonteCarlo Study
355	  const reco::GenParticle * muMc0 = mu0.genLepton();
356	  const reco::GenParticle * muMc1 = mu1.genLepton();
357	  const Candidate * motherMu0 =  muMc0->mother();
358	  const Candidate * motherMu1 =  muMc1->mother();
359	  int num_dau_muon0 = motherMu0->numberOfDaughters();
360	  int num_dau_muon1 = motherMu1->numberOfDaughters();
361	  cout<<"         muone0"<<endl;
362	  cout<<"         num di daughters = "<< num_dau_muon0 <<endl;
363	  if( num_dau_muon0 > 1 ){
364	    for(int j = 0; j <  num_dau_muon0; ++j){
365	      int id =motherMu0 ->daughter(j)->pdgId();
366	      cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
367	      if(id == 22) FSR_mu0=true;
368	    }
369	  }//end check of gamma
370
371	  cout<<"         muone1"<<endl;
372	  cout<<"         num di daughters = "<< num_dau_muon1 <<endl;
373	  if( num_dau_muon1 > 1 ){
374	    for(int j = 0; j <  num_dau_muon1; ++j){
375	      int id = motherMu1->daughter(j)->pdgId(); 
376	      cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
377	      if(id == 22) FSR_mu1=true;
378	    }
379	  }//end check of gamma
380	  if(FSR_mu0 || FSR_mu1 )h_zMuSamass_FSR->Fill(zmass);
381	  else h_zMuSamass_no_FSR->Fill(zmass);
382	  if(lep0->isGlobalMuon() && trig0found){ 
383	    if(FSR_mu1){
384	      h_staProbe_eta_FSR->Fill(eta1);
385	      h_staProbe_pt_FSR->Fill(pt1);
386	      h_staProbe_pt_eta_FSR_->Fill(pt1,eta1);
387
388	    }else{
389	      h_staProbe_eta_no_FSR->Fill(eta1);
390	      h_staProbe_pt_no_FSR->Fill(pt1);
391	      h_staProbe_pt_eta_no_FSR_->Fill(pt1,eta1);
392
393	    }
394	    if(FSR_mu1){
395	      h_Iso_FSR_->Fill(Tracker_isovalue_mu1);
396	      h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
397	    }
398	    else{ 
399	      h_Iso_->Fill(Tracker_isovalue_mu1);
400	      h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
401	    }
402	  }
403	  if(lep1->isGlobalMuon() && trig1found){ 
404	    if(FSR_mu0){
405	      h_staProbe_eta_FSR->Fill(eta0);
406	      h_staProbe_pt_FSR->Fill(pt0);
407	      h_staProbe_pt_eta_FSR_->Fill(pt0,eta0);
408
409	    }else{
410	      h_staProbe_eta_no_FSR->Fill(eta0);
411	      h_staProbe_pt_no_FSR->Fill(pt0);
412	      h_staProbe_pt_eta_FSR_->Fill(pt0,eta0);
413
414	    }
415	    if(FSR_mu0){
416	      h_Iso_FSR_->Fill(Tracker_isovalue_mu0);
417	      h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
418	    }
419	    else{ 
420	      h_Iso_->Fill(Tracker_isovalue_mu0);
421	      h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
422	    }
423	  }
424	}// end MC match
425      }//end of cuts
426     }// end loop on ZMuSa cand
427  }// end if ZMuSa size > 0
428
429  //ZMuTk  
430  if (zMuTk->size() > 0 ) {
431    event.getByLabel(zMuTkMatchMap_, zMuTkMatchMap); 
432    for(unsigned int i = 0; i < zMuTk->size(); ++i) { //loop on candidates
433      const Candidate & zMuTkCand = (*zMuTk)[i]; //the candidate
434      CandidateBaseRef zMuTkCandRef = zMuTk->refAt(i);
435      CandidateBaseRef dau0 = zMuTkCand.daughter(0)->masterClone();
436      CandidateBaseRef dau1 = zMuTkCand.daughter(1)->masterClone();
437      const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon
438      const pat::GenericParticle& mu1 = dynamic_cast<const pat::GenericParticle &>(*dau1);
439      
440      
441      double zmass= zMuTkCand.mass();
442      double pt0 = mu0.pt();
443      double pt1 = mu1.pt();
444      double eta0 = mu0.eta();
445      double eta1 = mu1.eta();
446      if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){//kinematical cuts
447	GenParticleRef zMuTkMatch = (*zMuTkMatchMap)[zMuTkCandRef];
448	if(zMuTkMatch.isNonnull()) {  // ZMuTk matched
449	  FSR_mu = false;
450	  FSR_tk = false;
451	  cout<<"          ZmuTk cuts && matched"<<endl;
452	  zmtcounter++;
453	  //Isodeposit
454	  const pat::IsoDeposit * muTrackIso =mu0.isoDeposit(pat::TrackIso);
455	  const pat::IsoDeposit * tkTrackIso =mu1.isoDeposit(pat::TrackIso);
456	  Direction muDir = Direction(mu0.eta(),mu0.phi());
457	  Direction tkDir = Direction(mu1.eta(),mu1.phi());
458	  
459	  IsoDeposit::AbsVetos vetos_mu;
460	  vetos_mu.push_back(new ConeVeto( muDir, dRVeto_ ));
461	  vetos_mu.push_back(new ThresholdVeto( ptThreshold_ ));
462	  
463	  reco::IsoDeposit::AbsVetos vetos_tk;
464	  vetos_tk.push_back(new ConeVeto( tkDir, dRVeto_ ));
465	  vetos_tk.push_back(new ThresholdVeto( ptThreshold_ ));
466	  
467	  double  Tracker_isovalue_mu = muTrackIso->sumWithin(dRTrk_,vetos_mu);
468	  double  Tracker_isovalue_tk = tkTrackIso->sumWithin(dRTrk_,vetos_tk);
469	  
470	  //MonteCarlo Study
471	  const reco::GenParticle * muMc0 = mu0.genLepton();
472	  const reco::GenParticle * muMc1 = mu1.genParticle() ;
473	  const Candidate * motherMu0 =  muMc0->mother();
474	  const Candidate * motherMu1 =  muMc1->mother();
475	  int num_dau_muon0 = motherMu0->numberOfDaughters();
476	  int num_dau_muon1 = motherMu1->numberOfDaughters();
477	  cout<<"numero di figli muone0 = " << num_dau_muon0 <<endl;
478	  cout<<"numero di figli muone1 = " << num_dau_muon1 <<endl; 
479	
480	  cout<<"         muon"<<endl;
481	  cout<<"         num di daughters = "<< num_dau_muon0 <<endl;
482	  if( num_dau_muon0 > 1 ){
483	    for(int j = 0; j <  num_dau_muon0; ++j){
484	      int id = motherMu0->daughter(j)->pdgId(); 
485	      cout<<"         dau["<<j<<"] pdg ID = "<<id<<endl;
486	      if(id == 22) {
487	 	FSR_mu=true;
488	      }
489	    }
490	  }//end check of gamma
491	  else cout<<"         dau[0] pdg ID = "<<motherMu0->daughter(0)->pdgId()<<endl;
492	  cout<<"         traccia"<<endl;
493	  cout<<"         num di daughters = "<< num_dau_muon1 <<endl;
494	  if( num_dau_muon1 > 1 ){
495	    for(int j = 0; j <  num_dau_muon1; ++j){
496	      int id = motherMu1->daughter(j)->pdgId(); 
497	      cout<<"         dau["<<j<<"] pdg ID = "<<id<<endl;
498	      if(id == 22) {
499		FSR_tk=true;
500	      }
501	    }
502	  }//end check of gamma
503	  else cout<<"         dau[0] pdg ID = "<<motherMu1->daughter(0)->pdgId()<<endl;
504	  cout<<"Mu Isolation = "<< Tracker_isovalue_mu <<endl;
505	  cout<<"Track Isolation = "<< Tracker_isovalue_tk <<endl;
506	  if(FSR_mu){
507	    h_Iso_FSR_->Fill(Tracker_isovalue_mu);
508	    h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu);
509	  }  
510	  else{ 
511	    h_Iso_->Fill( Tracker_isovalue_mu);
512	    h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu);
513
514	  }
515	  if(FSR_tk){
516	    h_Iso_FSR_->Fill(Tracker_isovalue_tk);
517	    h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_tk);
518	    h_trackProbe_eta_FSR->Fill(eta1);
519	    h_trackProbe_pt_FSR->Fill(pt1);
520 	  }
521	  else{
522	    h_Iso_->Fill( Tracker_isovalue_tk);
523	    h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_tk);
524	    h_trackProbe_eta_no_FSR->Fill(eta1);
525	    h_trackProbe_pt_no_FSR->Fill(pt1);
526	  }
527	}// end MC match
528      }//end Kine-cuts
529    }// end loop on ZMuTk cand
530  }// end if ZMuTk size > 0
531}// end analyze
532
533
534
535void ZMuMu_Radiative_analyzer::endJob() {
536  cout<<" ============= Summary =========="<<endl;
537  cout <<" Numero di eventi "<< evntcounter << endl; 
538  cout <<" 1)Numero di ZMuMu matched dopo i tagli cinematici = "<< zmmcounter << endl;
539  cout <<" 2)Numero di ZMuSa matched dopo i tagli cinematici = "<< zmscounter << endl;
540  cout <<" 3)Numero di ZMuTk matched dopo i tagli cinematici = "<< zmtcounter << endl;
541  double n1= h_Iso_FSR_->Integral();
542  double icut1= h_Iso_FSR_->Integral(0,15);
543  double eff_iso_FSR = (double)icut1/(double)n1;
544  double err_iso_FSR = sqrt(eff_iso_FSR*(1-eff_iso_FSR)/n1);
545  double n2= h_Iso_->Integral();
546  double icut2= h_Iso_->Integral(0,15);
547  double eff_iso= (double)icut2/(double)n2;
548  double err_iso = sqrt(eff_iso*(1-eff_iso)/n2);
549  cout<<" ============= Isolation Efficiecy =========="<<endl;
550  cout<<"Isolation Efficiency  = "<< eff_iso <<" +/- "<< err_iso <<endl;
551  cout<<"Isolation Efficiency with FSR = "<< eff_iso_FSR <<" +/- "<< err_iso_FSR <<endl;
552
553 }
554  
555#include "FWCore/Framework/interface/MakerMacros.h"
556
557DEFINE_FWK_MODULE(ZMuMu_Radiative_analyzer);
558