Reference Manual

Go to the documentation of this file.
00001 /* \class ZMuMu_Radiative_analyzer
00002  * 
00003  * author: Pasquale Noli
00004  *
00005  * ZMuMu Radiative analyzer
00006  * 
00007  *
00008  */
00009 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00010 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00011 #include "DataFormats/MuonReco/interface/Muon.h"
00012 #include "DataFormats/Common/interface/AssociationVector.h"
00013 #include "DataFormats/Candidate/interface/CandMatchMap.h" 
00014 #include "DataFormats/Candidate/interface/Particle.h"
00015 #include "DataFormats/Candidate/interface/Candidate.h"
00016 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00017 #include "DataFormats/Candidate/interface/OverlapChecker.h"
00018 #include "DataFormats/PatCandidates/interface/Muon.h" 
00019 #include "DataFormats/PatCandidates/interface/Lepton.h" 
00020 #include "DataFormats/PatCandidates/interface/GenericParticle.h"
00021 #include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h"
00022 #include "DataFormats/PatCandidates/interface/PATObject.h"
00023 #include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
00024 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
00025 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
00026 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
00027 #include "DataFormats/PatCandidates/interface/Isolation.h"
00028 #include "DataFormats/Common/interface/Handle.h"
00029 #include "DataFormats/Common/interface/ValueMap.h"
00030 #include "DataFormats/Candidate/interface/CandAssociation.h"
00031 #include "FWCore/Framework/interface/EDAnalyzer.h"
00032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00033 #include "FWCore/Framework/interface/Event.h"
00034 #include "FWCore/Framework/interface/EventSetup.h"
00035 #include "FWCore/Utilities/interface/InputTag.h"
00036 #include "FWCore/ServiceRegistry/interface/Service.h"
00037 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00038 #include <iostream>
00039 #include <iterator>
00040 #include <cmath>
00041 #include <vector>
00042 #include "TH1.h"
00043 #include "TH2.h"
00044 #include "TH3.h"
00045 
00046 
00047 using namespace edm;
00048 using namespace std;
00049 using namespace reco;
00050 using namespace isodeposit;
00051 
00052 class ZMuMu_Radiative_analyzer : public edm::EDAnalyzer {
00053 public:
00054   ZMuMu_Radiative_analyzer(const edm::ParameterSet& pset);
00055 private:
00056   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
00057   virtual void endJob();
00058 
00059   edm::InputTag zMuMu_, zMuMuMatchMap_, zMuTk_, zMuTkMatchMap_,zMuSa_, zMuSaMatchMap_;
00060   double dRVeto_, dRTrk_, ptThreshold_;
00061   //histograms 
00062   TH1D *h_zmass_FSR,*h_zmass_no_FSR;
00063   TH1D *h_zMuSamass_FSR,*h_zMuSamass_no_FSR;
00064   TH1D *h_zMuTkmass_FSR,*h_zMuTkmass_no_FSR;
00065   TH1D *h_Iso_,*h_Iso_FSR_ ;
00066   TH3D *h_Iso_3D_, *h_Iso_FSR_3D_; 
00067   TH2D *h_staProbe_pt_eta_no_FSR_, *h_staProbe_pt_eta_FSR_, *h_ProbeOk_pt_eta_no_FSR_, *h_ProbeOk_pt_eta_FSR_;
00068   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;
00069   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;
00070   //boolean 
00071   bool FSR_mu, FSR_tk, FSR_mu0, FSR_mu1;
00072   bool trig0found, trig1found;
00073   //counter
00074   int  zmmcounter , zmscounter, zmtcounter, evntcounter;
00075  };
00076 
00077 
00078 typedef edm::ValueMap<float> IsolationCollection;
00079 
00080 ZMuMu_Radiative_analyzer::ZMuMu_Radiative_analyzer(const ParameterSet& pset) : 
00081   zMuMu_(pset.getParameter<InputTag>("zMuMu")), 
00082   zMuMuMatchMap_(pset.getParameter<InputTag>("zMuMuMatchMap")),
00083   zMuTk_(pset.getParameter<InputTag>("zMuTk")), 
00084   zMuTkMatchMap_(pset.getParameter<InputTag>("zMuTkMatchMap")),
00085   zMuSa_(pset.getParameter<InputTag>("zMuSa")), 
00086   zMuSaMatchMap_(pset.getParameter<InputTag>("zMuSaMatchMap")),
00087   dRVeto_(pset.getUntrackedParameter<double>("veto")),
00088   dRTrk_(pset.getUntrackedParameter<double>("deltaRTrk")),
00089   ptThreshold_(pset.getUntrackedParameter<double>("ptThreshold")){ 
00090   zmmcounter=0;
00091   zmscounter=0;
00092   zmtcounter=0;
00093   evntcounter=0;
00094   Service<TFileService> fs;
00095    
00096   // general histograms
00097   h_zmass_FSR= fs->make<TH1D>("h_zmass_FRS","Invariant Z mass distribution",200,0,200);
00098   h_zmass_no_FSR= fs->make<TH1D>("h_zmass_no_FSR","Invariant Z mass distribution",200,0,200);
00099   h_zMuSamass_FSR= fs->make<TH1D>("h_zMuSamass_FRS","Invariant Z mass distribution",200,0,200);
00100   h_zMuSamass_no_FSR= fs->make<TH1D>("h_zMuSamass_no_FSR","Invariant Z mass distribution",200,0,200);
00101   h_zMuTkmass_FSR= fs->make<TH1D>("h_zMuTkmass_FRS","Invariant Z mass distribution",200,0,200);
00102   h_zMuTkmass_no_FSR= fs->make<TH1D>("h_zMuTkmass_no_FSR","Invariant Z mass distribution",200,0,200);
00103   h_Iso_= fs->make<TH1D>("h_iso","Isolation distribution of muons without FSR",100,0,20);
00104   h_Iso_FSR_= fs->make<TH1D>("h_iso_FSR","Isolation distribution of muons with FSR ",100,0,20);
00105   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);
00106   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);
00107   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);
00108   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);
00109   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);
00110   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);
00111   
00112   h_trackProbe_eta_no_FSR = fs->make<TH1D>("trackProbeEta_no_FSR","Eta of tracks",100,-2.0,2.0);
00113   h_trackProbe_pt_no_FSR = fs->make<TH1D>("trackProbePt_no_FSR","Pt of tracks",100,0,100);
00114   h_staProbe_eta_no_FSR = fs->make<TH1D>("standAloneProbeEta_no_FSR","Eta of standAlone",100,-2.0,2.0);
00115   h_staProbe_pt_no_FSR = fs->make<TH1D>("standAloneProbePt_no_FSR","Pt of standAlone",100,0,100);
00116   h_ProbeOk_eta_no_FSR = fs->make<TH1D>("probeOkEta_no_FSR","Eta of probe Ok",100,-2.0,2.0);
00117   h_ProbeOk_pt_no_FSR = fs->make<TH1D>("probeOkPt_no_FSR","Pt of probe ok",100,0,100);
00118 
00119   h_trackProbe_eta_FSR = fs->make<TH1D>("trackProbeEta_FSR","Eta of tracks",100,-2.0,2.0);
00120   h_trackProbe_pt_FSR  = fs->make<TH1D>("trackProbePt_FSR","Pt of tracks",100,0,100);
00121   h_staProbe_eta_FSR  = fs->make<TH1D>("standAloneProbeEta_FSR","Eta of standAlone",100,-2.0,2.0);
00122   h_staProbe_pt_FSR  = fs->make<TH1D>("standAloneProbePt_FSR","Pt of standAlone",100,0,100);
00123   h_ProbeOk_eta_FSR  = fs->make<TH1D>("probeOkEta_FSR","Eta of probe Ok",100,-2.0,2.0);
00124   h_ProbeOk_pt_FSR  = fs->make<TH1D>("probeOkPt_FSR","Pt of probe ok",100,0,100);
00125 
00126 
00127 
00128 }
00129 
00130 void ZMuMu_Radiative_analyzer::analyze(const Event& event, const EventSetup& setup) {
00131   evntcounter++;
00132   Handle<CandidateView> zMuMu;                //Collection of Z made by  Mu global + Mu global 
00133   Handle<GenParticleMatch> zMuMuMatchMap;     //Map of Z made by Mu global + Mu global with MC
00134   event.getByLabel(zMuMu_, zMuMu); 
00135   Handle<CandidateView> zMuTk;                //Collection of Z made by  Mu global + Track 
00136   Handle<GenParticleMatch> zMuTkMatchMap;
00137   event.getByLabel(zMuTk_, zMuTk); 
00138   Handle<CandidateView> zMuSa;                //Collection of Z made by  Mu global + Sa 
00139   Handle<GenParticleMatch> zMuSaMatchMap;
00140   event.getByLabel(zMuSa_, zMuSa); 
00141   cout << "**********  New Event  ***********"<<endl; 
00142   // ZMuMu
00143   if (zMuMu->size() > 0 ) {
00144     event.getByLabel(zMuMuMatchMap_, zMuMuMatchMap); 
00145      for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates
00146      
00147       const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate
00148       CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i);
00149 
00150        
00151       CandidateBaseRef dau0 = zMuMuCand.daughter(0)->masterClone();
00152       CandidateBaseRef dau1 = zMuMuCand.daughter(1)->masterClone();
00153       const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon
00154       const pat::Muon& mu1 = dynamic_cast<const pat::Muon&>(*dau1);
00155             
00156       double zmass= zMuMuCand.mass();
00157       double pt0 = mu0.pt();
00158       double pt1 = mu1.pt();
00159       double eta0 = mu0.eta();
00160       double eta1 = mu1.eta();
00161       if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){
00162         GenParticleRef zMuMuMatch = (*zMuMuMatchMap)[zMuMuCandRef];
00163         if(zMuMuMatch.isNonnull()) {  // ZMuMu matched
00164           zmmcounter++;
00165           cout<<"         Zmumu cuts && matched" <<endl;
00166           FSR_mu0 = false;
00167           FSR_mu1 = false;
00168           
00169           //Isodeposit
00170           const pat::IsoDeposit * mu0TrackIso =mu0.isoDeposit(pat::TrackIso);
00171           const pat::IsoDeposit * mu1TrackIso =mu1.isoDeposit(pat::TrackIso);
00172           Direction mu0Dir = Direction(mu0.eta(),mu0.phi());
00173           Direction mu1Dir = Direction(mu1.eta(),mu1.phi());
00174           
00175           reco::IsoDeposit::AbsVetos vetos_mu0;
00176           vetos_mu0.push_back(new ConeVeto( mu0Dir, dRVeto_ ));
00177           vetos_mu0.push_back(new ThresholdVeto( ptThreshold_ ));
00178           
00179           reco::IsoDeposit::AbsVetos vetos_mu1;
00180           vetos_mu1.push_back(new ConeVeto( mu1Dir, dRVeto_ ));
00181           vetos_mu1.push_back(new ThresholdVeto( ptThreshold_ ));
00182           
00183           double  Tracker_isovalue_mu0 = mu0TrackIso->sumWithin(dRTrk_,vetos_mu0);
00184           double  Tracker_isovalue_mu1 = mu1TrackIso->sumWithin(dRTrk_,vetos_mu1);
00185          
00186           //trigger study 
00187           trig0found = false;
00188           trig1found = false;
00189           const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
00190             mu0.triggerObjectMatchesByPath( "HLT_Mu9" );
00191           const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 
00192             mu1.triggerObjectMatchesByPath( "HLT_Mu9" );
00193           if( mu0HLTMatches.size()>0 )
00194             trig0found = true;
00195           if( mu1HLTMatches.size()>0 )
00196             trig1found = true;
00197       
00198           //MonteCarlo Study
00199           const reco::GenParticle * muMc0 = mu0.genLepton();
00200           const reco::GenParticle * muMc1 = mu1.genLepton();
00201           const Candidate * motherMu0 =  muMc0->mother();
00202           int num_dau_muon0 = motherMu0->numberOfDaughters();
00203           const Candidate * motherMu1 =  muMc1->mother();
00204           int num_dau_muon1 = motherMu1->numberOfDaughters();
00205           cout<<"         muone0"<<endl;
00206           cout<<"         num di daughters = "<< num_dau_muon0 <<endl;
00207           if( num_dau_muon0 > 1 ){
00208             for(int j = 0; j <  num_dau_muon0; ++j){
00209               int id =motherMu0 ->daughter(j)->pdgId();
00210               cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
00211               if(id == 22) FSR_mu0=true;
00212             }
00213           }//end check of gamma
00214 
00215           cout<<"         muone1"<<endl;
00216           cout<<"         num di daughters = "<< num_dau_muon1 <<endl;
00217           if( num_dau_muon1 > 1 ){
00218             for(int j = 0; j <  num_dau_muon1; ++j){
00219               int id = motherMu1->daughter(j)->pdgId(); 
00220               cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
00221               if(id == 22) FSR_mu1=true;
00222             }
00223           }//end check of gamma
00224         
00225           if(FSR_mu0 || FSR_mu1 )h_zmass_FSR->Fill(zmass);
00226           else h_zmass_no_FSR->Fill(zmass);
00227 
00228           if (trig1found) {       // check efficiency of muon0 not imposing the trigger on it 
00229             cout<<"muon 1 is triggered "<<endl;
00230             if(FSR_mu0){
00231               cout<< "and muon 0 does FSR"<<endl;
00232               h_trackProbe_eta_FSR->Fill(eta0);
00233               h_trackProbe_pt_FSR->Fill(pt0);
00234               h_staProbe_eta_FSR->Fill(eta0);
00235               h_staProbe_pt_FSR->Fill(pt0);
00236               h_staProbe_pt_eta_FSR_->Fill(pt0,eta0);
00237               h_ProbeOk_eta_FSR->Fill(eta0);
00238               h_ProbeOk_pt_FSR->Fill(pt0);
00239               h_ProbeOk_pt_eta_FSR_->Fill(pt0,eta0);
00240             }else{
00241               cout<<"and muon 0 doesn't FSR"<<endl;
00242               h_trackProbe_eta_no_FSR->Fill(eta0);
00243               h_trackProbe_pt_no_FSR->Fill(pt0);
00244               h_staProbe_eta_no_FSR->Fill(eta0);
00245               h_staProbe_pt_no_FSR->Fill(pt0);
00246               h_staProbe_pt_eta_no_FSR_->Fill(pt0,eta0);
00247               h_ProbeOk_eta_no_FSR->Fill(eta0);
00248               h_ProbeOk_pt_no_FSR->Fill(pt0);
00249               h_ProbeOk_pt_eta_no_FSR_->Fill(pt0,eta0);
00250             }
00251             if(FSR_mu0){
00252               h_Iso_FSR_->Fill(Tracker_isovalue_mu0);
00253               h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
00254             }
00255             else{
00256               h_Iso_->Fill(Tracker_isovalue_mu0);
00257               h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
00258             }
00259           }
00260           if (trig0found) {         // check efficiency of muon1 not imposing the trigger on it 
00261             cout<<"muon 0 is triggered"<<endl;
00262             if(FSR_mu1){ 
00263               cout<<"and muon 1  does FSR"<<endl;
00264               h_trackProbe_eta_FSR->Fill(eta1);
00265               h_staProbe_eta_FSR->Fill(eta1);
00266               h_trackProbe_pt_FSR->Fill(pt1);
00267               h_staProbe_pt_FSR->Fill(pt1);
00268               h_ProbeOk_eta_FSR->Fill(eta1);
00269               h_ProbeOk_pt_FSR->Fill(pt1);
00270               h_staProbe_pt_eta_FSR_->Fill(pt1,eta1);
00271               h_ProbeOk_pt_eta_FSR_->Fill(pt1,eta1);
00272 
00273             }else{
00274               cout<<"and muon 1 doesn't FSR"<<endl;
00275               h_trackProbe_eta_no_FSR->Fill(eta1);
00276               h_staProbe_eta_no_FSR->Fill(eta1);
00277               h_trackProbe_pt_no_FSR->Fill(pt1);
00278               h_staProbe_pt_no_FSR->Fill(pt1);
00279               h_ProbeOk_eta_no_FSR->Fill(eta1);
00280               h_ProbeOk_pt_no_FSR->Fill(pt1);
00281               h_staProbe_pt_eta_no_FSR_->Fill(pt1,eta1);
00282               h_ProbeOk_pt_eta_no_FSR_->Fill(pt1,eta1);
00283 
00284 
00285             }
00286             if(FSR_mu1){
00287               h_Iso_FSR_->Fill(Tracker_isovalue_mu1);
00288               h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
00289             }else{ 
00290               h_Iso_->Fill(Tracker_isovalue_mu1);
00291               h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
00292             }
00293           }
00294         }// end MC match
00295       }//end of cuts
00296      }// end loop on ZMuMu cand
00297   }// end if ZMuMu size > 0
00298   
00299   // ZMuSa
00300   if (zMuSa->size() > 0 ) {
00301     event.getByLabel(zMuSaMatchMap_, zMuSaMatchMap); 
00302      for(unsigned int i = 0; i < zMuSa->size(); ++i) { //loop on candidates
00303      
00304       const Candidate & zMuSaCand = (*zMuSa)[i]; //the candidate
00305       CandidateBaseRef zMuSaCandRef = zMuSa->refAt(i);
00306       const Candidate *  lep0 =zMuSaCand.daughter(0);  
00307       const Candidate *  lep1 =zMuSaCand.daughter(1);  
00308       CandidateBaseRef dau0 = lep0->masterClone();
00309       CandidateBaseRef dau1 = lep1->masterClone();
00310       const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon
00311       const pat::Muon& mu1 = dynamic_cast<const pat::Muon&>(*dau1);
00312             
00313       double zmass= zMuSaCand.mass();
00314       double pt0 = mu0.pt();
00315       double pt1 = mu1.pt();
00316       double eta0 = mu0.eta();
00317       double eta1 = mu1.eta();
00318       if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){
00319         GenParticleRef zMuSaMatch = (*zMuSaMatchMap)[zMuSaCandRef];
00320         if(zMuSaMatch.isNonnull()) {  // ZMuSa matched
00321           cout<<"         Zmusa cuts && matched" <<endl;
00322           FSR_mu0 = false;
00323           FSR_mu1 = false;
00324           zmscounter++;   
00325           //Isodeposit
00326           const pat::IsoDeposit * mu0TrackIso =mu0.isoDeposit(pat::TrackIso);
00327           const pat::IsoDeposit * mu1TrackIso =mu1.isoDeposit(pat::TrackIso);
00328           Direction mu0Dir = Direction(mu0.eta(),mu0.phi());
00329           Direction mu1Dir = Direction(mu1.eta(),mu1.phi());
00330           
00331           reco::IsoDeposit::AbsVetos vetos_mu0;
00332           vetos_mu0.push_back(new ConeVeto( mu0Dir, dRVeto_ ));
00333           vetos_mu0.push_back(new ThresholdVeto( ptThreshold_ ));
00334           
00335           reco::IsoDeposit::AbsVetos vetos_mu1;
00336           vetos_mu1.push_back(new ConeVeto( mu1Dir, dRVeto_ ));
00337           vetos_mu1.push_back(new ThresholdVeto( ptThreshold_ ));
00338           
00339           double  Tracker_isovalue_mu0 = mu0TrackIso->sumWithin(dRTrk_,vetos_mu0);
00340           double  Tracker_isovalue_mu1 = mu1TrackIso->sumWithin(dRTrk_,vetos_mu1);
00341 
00342           // HLT match (check just dau0 the global)
00343           const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 
00344             mu0.triggerObjectMatchesByPath( "HLT_Mu9" );
00345           const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 
00346             mu1.triggerObjectMatchesByPath( "HLT_Mu9" );
00347           trig0found = false;
00348           trig1found = false;
00349           if( mu0HLTMatches.size()>0 )
00350             trig0found = true;
00351           if( mu1HLTMatches.size()>0 )
00352             trig1found = true;
00353           
00354           //MonteCarlo Study
00355           const reco::GenParticle * muMc0 = mu0.genLepton();
00356           const reco::GenParticle * muMc1 = mu1.genLepton();
00357           const Candidate * motherMu0 =  muMc0->mother();
00358           const Candidate * motherMu1 =  muMc1->mother();
00359           int num_dau_muon0 = motherMu0->numberOfDaughters();
00360           int num_dau_muon1 = motherMu1->numberOfDaughters();
00361           cout<<"         muone0"<<endl;
00362           cout<<"         num di daughters = "<< num_dau_muon0 <<endl;
00363           if( num_dau_muon0 > 1 ){
00364             for(int j = 0; j <  num_dau_muon0; ++j){
00365               int id =motherMu0 ->daughter(j)->pdgId();
00366               cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
00367               if(id == 22) FSR_mu0=true;
00368             }
00369           }//end check of gamma
00370 
00371           cout<<"         muone1"<<endl;
00372           cout<<"         num di daughters = "<< num_dau_muon1 <<endl;
00373           if( num_dau_muon1 > 1 ){
00374             for(int j = 0; j <  num_dau_muon1; ++j){
00375               int id = motherMu1->daughter(j)->pdgId(); 
00376               cout<<"         dauther["<<j<<"] pdgId = "<<id<<endl; 
00377               if(id == 22) FSR_mu1=true;
00378             }
00379           }//end check of gamma
00380           if(FSR_mu0 || FSR_mu1 )h_zMuSamass_FSR->Fill(zmass);
00381           else h_zMuSamass_no_FSR->Fill(zmass);
00382           if(lep0->isGlobalMuon() && trig0found){ 
00383             if(FSR_mu1){
00384               h_staProbe_eta_FSR->Fill(eta1);
00385               h_staProbe_pt_FSR->Fill(pt1);
00386               h_staProbe_pt_eta_FSR_->Fill(pt1,eta1);
00387 
00388             }else{
00389               h_staProbe_eta_no_FSR->Fill(eta1);
00390               h_staProbe_pt_no_FSR->Fill(pt1);
00391               h_staProbe_pt_eta_no_FSR_->Fill(pt1,eta1);
00392 
00393             }
00394             if(FSR_mu1){
00395               h_Iso_FSR_->Fill(Tracker_isovalue_mu1);
00396               h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
00397             }
00398             else{ 
00399               h_Iso_->Fill(Tracker_isovalue_mu1);
00400               h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1);
00401             }
00402           }
00403           if(lep1->isGlobalMuon() && trig1found){ 
00404             if(FSR_mu0){
00405               h_staProbe_eta_FSR->Fill(eta0);
00406               h_staProbe_pt_FSR->Fill(pt0);
00407               h_staProbe_pt_eta_FSR_->Fill(pt0,eta0);
00408 
00409             }else{
00410               h_staProbe_eta_no_FSR->Fill(eta0);
00411               h_staProbe_pt_no_FSR->Fill(pt0);
00412               h_staProbe_pt_eta_FSR_->Fill(pt0,eta0);
00413 
00414             }
00415             if(FSR_mu0){
00416               h_Iso_FSR_->Fill(Tracker_isovalue_mu0);
00417               h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
00418             }
00419             else{ 
00420               h_Iso_->Fill(Tracker_isovalue_mu0);
00421               h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0);
00422             }
00423           }
00424         }// end MC match
00425       }//end of cuts
00426      }// end loop on ZMuSa cand
00427   }// end if ZMuSa size > 0
00428 
00429   //ZMuTk  
00430   if (zMuTk->size() > 0 ) {
00431     event.getByLabel(zMuTkMatchMap_, zMuTkMatchMap); 
00432     for(unsigned int i = 0; i < zMuTk->size(); ++i) { //loop on candidates
00433       const Candidate & zMuTkCand = (*zMuTk)[i]; //the candidate
00434       CandidateBaseRef zMuTkCandRef = zMuTk->refAt(i);
00435       CandidateBaseRef dau0 = zMuTkCand.daughter(0)->masterClone();
00436       CandidateBaseRef dau1 = zMuTkCand.daughter(1)->masterClone();
00437       const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon
00438       const pat::GenericParticle& mu1 = dynamic_cast<const pat::GenericParticle &>(*dau1);
00439       
00440       
00441       double zmass= zMuTkCand.mass();
00442       double pt0 = mu0.pt();
00443       double pt1 = mu1.pt();
00444       double eta0 = mu0.eta();
00445       double eta1 = mu1.eta();
00446       if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){//kinematical cuts
00447         GenParticleRef zMuTkMatch = (*zMuTkMatchMap)[zMuTkCandRef];
00448         if(zMuTkMatch.isNonnull()) {  // ZMuTk matched
00449           FSR_mu = false;
00450           FSR_tk = false;
00451           cout<<"          ZmuTk cuts && matched"<<endl;
00452           zmtcounter++;
00453           //Isodeposit
00454           const pat::IsoDeposit * muTrackIso =mu0.isoDeposit(pat::TrackIso);
00455           const pat::IsoDeposit * tkTrackIso =mu1.isoDeposit(pat::TrackIso);
00456           Direction muDir = Direction(mu0.eta(),mu0.phi());
00457           Direction tkDir = Direction(mu1.eta(),mu1.phi());
00458           
00459           IsoDeposit::AbsVetos vetos_mu;
00460           vetos_mu.push_back(new ConeVeto( muDir, dRVeto_ ));
00461           vetos_mu.push_back(new ThresholdVeto( ptThreshold_ ));
00462           
00463           reco::IsoDeposit::AbsVetos vetos_tk;
00464           vetos_tk.push_back(new ConeVeto( tkDir, dRVeto_ ));
00465           vetos_tk.push_back(new ThresholdVeto( ptThreshold_ ));
00466           
00467           double  Tracker_isovalue_mu = muTrackIso->sumWithin(dRTrk_,vetos_mu);
00468           double  Tracker_isovalue_tk = tkTrackIso->sumWithin(dRTrk_,vetos_tk);
00469           
00470           //MonteCarlo Study
00471           const reco::GenParticle * muMc0 = mu0.genLepton();
00472           const reco::GenParticle * muMc1 = mu1.genParticle() ;
00473           const Candidate * motherMu0 =  muMc0->mother();
00474           const Candidate * motherMu1 =  muMc1->mother();
00475           int num_dau_muon0 = motherMu0->numberOfDaughters();
00476           int num_dau_muon1 = motherMu1->numberOfDaughters();
00477           cout<<"numero di figli muone0 = " << num_dau_muon0 <<endl;
00478           cout<<"numero di figli muone1 = " << num_dau_muon1 <<endl; 
00479         
00480           cout<<"         muon"<<endl;
00481           cout<<"         num di daughters = "<< num_dau_muon0 <<endl;
00482           if( num_dau_muon0 > 1 ){
00483             for(int j = 0; j <  num_dau_muon0; ++j){
00484               int id = motherMu0->daughter(j)->pdgId(); 
00485               cout<<"         dau["<<j<<"] pdg ID = "<<id<<endl;
00486               if(id == 22) {
00487                 FSR_mu=true;
00488               }
00489             }
00490           }//end check of gamma
00491           else cout<<"         dau[0] pdg ID = "<<motherMu0->daughter(0)->pdgId()<<endl;
00492           cout<<"         traccia"<<endl;
00493           cout<<"         num di daughters = "<< num_dau_muon1 <<endl;
00494           if( num_dau_muon1 > 1 ){
00495             for(int j = 0; j <  num_dau_muon1; ++j){
00496               int id = motherMu1->daughter(j)->pdgId(); 
00497               cout<<"         dau["<<j<<"] pdg ID = "<<id<<endl;
00498               if(id == 22) {
00499                 FSR_tk=true;
00500               }
00501             }
00502           }//end check of gamma
00503           else cout<<"         dau[0] pdg ID = "<<motherMu1->daughter(0)->pdgId()<<endl;
00504           cout<<"Mu Isolation = "<< Tracker_isovalue_mu <<endl;
00505           cout<<"Track Isolation = "<< Tracker_isovalue_tk <<endl;
00506           if(FSR_mu){
00507             h_Iso_FSR_->Fill(Tracker_isovalue_mu);
00508             h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu);
00509           }  
00510           else{ 
00511             h_Iso_->Fill( Tracker_isovalue_mu);
00512             h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu);
00513 
00514           }
00515           if(FSR_tk){
00516             h_Iso_FSR_->Fill(Tracker_isovalue_tk);
00517             h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_tk);
00518             h_trackProbe_eta_FSR->Fill(eta1);
00519             h_trackProbe_pt_FSR->Fill(pt1);
00520           }
00521           else{
00522             h_Iso_->Fill( Tracker_isovalue_tk);
00523             h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_tk);
00524             h_trackProbe_eta_no_FSR->Fill(eta1);
00525             h_trackProbe_pt_no_FSR->Fill(pt1);
00526           }
00527         }// end MC match
00528       }//end Kine-cuts
00529     }// end loop on ZMuTk cand
00530   }// end if ZMuTk size > 0
00531 }// end analyze
00532 
00533 
00534 
00535 void ZMuMu_Radiative_analyzer::endJob() {
00536   cout<<" ============= Summary =========="<<endl;
00537   cout <<" Numero di eventi "<< evntcounter << endl; 
00538   cout <<" 1)Numero di ZMuMu matched dopo i tagli cinematici = "<< zmmcounter << endl;
00539   cout <<" 2)Numero di ZMuSa matched dopo i tagli cinematici = "<< zmscounter << endl;
00540   cout <<" 3)Numero di ZMuTk matched dopo i tagli cinematici = "<< zmtcounter << endl;
00541   double n1= h_Iso_FSR_->Integral();
00542   double icut1= h_Iso_FSR_->Integral(0,15);
00543   double eff_iso_FSR = (double)icut1/(double)n1;
00544   double err_iso_FSR = sqrt(eff_iso_FSR*(1-eff_iso_FSR)/n1);
00545   double n2= h_Iso_->Integral();
00546   double icut2= h_Iso_->Integral(0,15);
00547   double eff_iso= (double)icut2/(double)n2;
00548   double err_iso = sqrt(eff_iso*(1-eff_iso)/n2);
00549   cout<<" ============= Isolation Efficiecy =========="<<endl;
00550   cout<<"Isolation Efficiency  = "<< eff_iso <<" +/- "<< err_iso <<endl;
00551   cout<<"Isolation Efficiency with FSR = "<< eff_iso_FSR <<" +/- "<< err_iso_FSR <<endl;
00552 
00553  }
00554   
00555 #include "FWCore/Framework/interface/MakerMacros.h"
00556 
00557 DEFINE_FWK_MODULE(ZMuMu_Radiative_analyzer);
00558
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMu_radiative_analysis.cc
