Reference Manual

Go to the documentation of this file.
00001 #include "DataFormats/Common/interface/AssociationVector.h"
00002 #include "FWCore/Framework/interface/EDAnalyzer.h"
00003 #include "DataFormats/Candidate/interface/Candidate.h"
00004 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00006 #include "FWCore/Framework/interface/Event.h"
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 #include "FWCore/Utilities/interface/InputTag.h"
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00011 #include "DataFormats/Math/interface/deltaR.h"
00012 #include "DataFormats/PatCandidates/interface/Muon.h"
00013 #include "DataFormats/TrackReco/interface/Track.h"
00014 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00015 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
00016 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
00017 #include "DataFormats/PatCandidates/interface/Isolation.h"
00018 #include "DataFormats/Common/interface/ValueMap.h"
00019 #include "DataFormats/PatCandidates/interface/GenericParticle.h"
00020 #include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
00021 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
00022 #include "TH1.h"
00023 #include "TH2.h"
00024 #include "TMath.h"
00025 #include <vector>
00026 #include <string>
00027 #include <iostream>
00028 #include <sstream>
00029 
00030 using namespace edm;
00031 using namespace std;
00032 using namespace reco;
00033 using namespace isodeposit;
00034 
00035 class ZMuMuIsolationAnalyzer : public edm::EDAnalyzer {
00036 public:
00037   ZMuMuIsolationAnalyzer(const edm::ParameterSet& pset);
00038 private:
00039   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
00040   virtual void endJob();
00041   InputTag src_muons;
00042   double dRVeto;
00043   double dRTrk, dREcal, dRHcal;
00044   double ptThreshold, etEcalThreshold, etHcalThreshold;
00045   double alpha, beta;
00046   double pt, eta;
00047   double iso_cut;
00048 
00049   TH1F * h_IsoZ_tk,* h_IsoW_tk,* h_IsoOther_tk ;
00050   TH1F * h_IsoZ_ecal,* h_IsoW_ecal,* h_IsoOther_ecal ;
00051   TH1F * h_IsoZ_hcal,* h_IsoW_hcal,* h_IsoOther_hcal ;
00052   TH1F * IsoZ,* IsoW,* IsoOther ;
00053   TH1F * TkrPt,* EcalEt,* HcalEt ;
00054   TH1F * EcalEtZ, * HcalEtZ;
00055 
00056   TH1F * Z_eta,* W_eta,* Other_eta;
00057   TH1F * Z_eta_postSelection,* W_eta_postSelection,* Other_eta_postSelection;
00058   TH1F * Z_pt,* W_pt,* Other_pt;
00059   TH1F * Z_pt_postSelection,* W_pt_postSelection,* Other_pt_postSelection;
00060 
00061   enum MuTag { muFromZ, muFromW, muFromOther };
00062   template<typename T>
00063   MuTag muTag(const T & mu) const;
00064   void Deposits(const pat::IsoDeposit* isodep, double dR_max,  TH1F* hist);
00065   void histo(const TH1F* hist, const char* cx, const char* cy) const;
00066 };
00067 
00068 template<typename T>
00069 ZMuMuIsolationAnalyzer::MuTag ZMuMuIsolationAnalyzer::muTag(const T& mu) const {
00070   GenParticleRef p = mu.genParticleRef();
00071   if(p.isNull()){
00072     //   cout<<"genParticleRef is null "<<endl;
00073     return muFromOther;
00074 }
00075   int sizem = p->numberOfMothers();
00076   if(sizem != 1) {
00077     //cout<<"number of mothers !=1 "<<endl;
00078  return muFromOther;
00079   }
00080   const Candidate * moth1 = p->mother();
00081   if(moth1 == 0) {
00082     return muFromOther;
00083     //cout<<"no mother "<<endl;
00084   }
00085   int pdgId1 = moth1->pdgId();
00086   if(abs(pdgId1)!=13){ 
00087     return muFromOther;
00088     //cout<<"mother is not a muon"<<endl;
00089   }
00090   const Candidate * moth2 = moth1->mother();
00091   if(moth2 == 0) {
00092     return muFromOther;
00093     //cout<<"no mother "<<endl;
00094 }
00095   int pdgId2 = moth2->pdgId();
00096   if(pdgId2 == 23) {
00097     //cout<<" muon from Z"<<endl;
00098     return muFromZ;
00099   }
00100   if(abs(pdgId2)==24) return muFromW;
00101   else {
00102     //cout<<" muon from other"<<endl;
00103     return muFromOther;
00104   }
00105 }
00106 
00107 void ZMuMuIsolationAnalyzer::Deposits(const pat::IsoDeposit* isodep,double dR_max,TH1F* hist){
00108   for(IsoDeposit::const_iterator it= isodep->begin(); it!= isodep->end(); ++it){
00109     if(it->dR()<dR_max) {
00110       double theta= 2*(TMath::ATan(TMath::Exp(-(it->eta() ) ) ) );
00111       // double theta= 2;
00112       hist->Fill(it->value()/TMath::Sin(theta));
00113     }
00114   }
00115 }
00116 
00117 void ZMuMuIsolationAnalyzer::histo(const TH1F* hist, const char* cx, const char*cy) const{
00118   hist->GetXaxis()->SetTitle(cx);
00119   hist->GetYaxis()->SetTitle(cy);
00120   hist->GetXaxis()->SetTitleOffset(1);
00121   hist->GetYaxis()->SetTitleOffset(1.2);
00122   hist->GetXaxis()->SetTitleSize(0.04);
00123   hist->GetYaxis()->SetTitleSize(0.04);
00124   hist->GetXaxis()->SetLabelSize(0.03);
00125   hist->GetYaxis()->SetLabelSize(0.03);
00126 }
00127 
00128 ZMuMuIsolationAnalyzer::ZMuMuIsolationAnalyzer(const ParameterSet& pset):
00129   src_muons(pset.getParameter<InputTag>("src")),
00130   dRVeto(pset.getUntrackedParameter<double>("veto")),
00131   dRTrk(pset.getUntrackedParameter<double>("deltaRTrk")),
00132   dREcal(pset.getUntrackedParameter<double>("deltaREcal")),
00133   dRHcal(pset.getUntrackedParameter<double>("deltaRHcal")),
00134   ptThreshold(pset.getUntrackedParameter<double>("ptThreshold")),
00135   etEcalThreshold(pset.getUntrackedParameter<double>("etEcalThreshold")),
00136   etHcalThreshold(pset.getUntrackedParameter<double>("etHcalThreshold")),
00137   alpha(pset.getUntrackedParameter<double>("alpha")),
00138   beta(pset.getUntrackedParameter<double>("beta")),
00139   pt(pset.getUntrackedParameter<double>("pt")),
00140   eta(pset.getUntrackedParameter<double>("eta")),
00141   iso_cut(pset.getUntrackedParameter<double>("isoCut")) {
00142   edm::Service<TFileService> fs;
00143   std::ostringstream str1,str2,str3,str4,str5,str6,str7,str8,str9,str10,n_tracks;
00144   str1 << "muons from Z with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk;
00145   str2 << "muons from W with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk;
00146   str3 << "muons from Others with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk;
00147   str4 << "muons from Z with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
00148   str5 << "muons from W with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
00149   str6 << "muons from Other with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
00150   n_tracks <<"Number of tracks for muon with p_{t} > " << ptThreshold <<" and #Delta R < "<<dRTrk<< " GeV/c";
00151   str7<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<"(Tracker)";
00152   str8<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<"(Ecal)";
00153   str9<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<"(Hcal)";
00154   str10<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
00155   h_IsoZ_tk = fs->make<TH1F>("ZIso_Tk",str1.str().c_str(),100,0.,20.);
00156   h_IsoW_tk = fs->make<TH1F>("WIso_Tk",str2.str().c_str(),100,0.,20.);
00157   h_IsoOther_tk = fs->make<TH1F>("otherIso_Tk",str3.str().c_str(),100,0.,20.);
00158   h_IsoZ_ecal = fs->make<TH1F>("ZIso_ecal",str1.str().c_str(),100,0.,20.);
00159   h_IsoW_ecal = fs->make<TH1F>("WIso_ecal",str2.str().c_str(),100,0.,20.);
00160   h_IsoOther_ecal = fs->make<TH1F>("otherIso_ecal",str3.str().c_str(),100,0.,20.);
00161   h_IsoZ_hcal = fs->make<TH1F>("ZIso_hcal",str1.str().c_str(),100,0.,20.);
00162   h_IsoW_hcal = fs->make<TH1F>("WIso_hcal",str2.str().c_str(),100,0.,20.);
00163   h_IsoOther_hcal = fs->make<TH1F>("otherIso_hcal",str3.str().c_str(),100,0.,20.);
00164   IsoZ = fs->make<TH1F>("ZIso",str4.str().c_str(),100,0.,20.);
00165   IsoW = fs->make<TH1F>("WIso",str5.str().c_str(),100,0.,20.);
00166   IsoOther = fs->make<TH1F>("otherIso",str6.str().c_str(),100,0.,20.);
00167 
00168  
00169   Z_eta = fs->make<TH1F>("Z_eta","#eta distribution for muons coming from Z",40,-eta,eta);
00170   W_eta = fs->make<TH1F>("W_eta","#eta distribution for muons coming from W",40,-eta,eta);
00171   Other_eta = fs->make<TH1F>("Other_eta","#eta distribution for muons coming from other",40,-eta,eta);
00172   Z_eta_postSelection = fs->make<TH1F>("Z_eta_postSelection","#eta distribution for muons coming from Z after iso selection",40,-eta,eta);
00173   W_eta_postSelection = fs->make<TH1F>("W_eta_postSelection","#eta distribution for muons coming from W after iso selection",40,-eta,eta);
00174   Other_eta_postSelection = fs->make<TH1F>("Other_eta_postSelection","#eta distribution for muons coming from other after iso selection",40,-eta,eta);
00175   
00176   Z_pt = fs->make<TH1F>("Z_pt","p_{T} distribution for muons coming from Z",40,pt,150.);
00177   W_pt = fs->make<TH1F>("W_pt","p_{T} distribution for muons coming from W",40,pt,150.);
00178   Other_pt = fs->make<TH1F>("Other_pt","p_{T} distribution for muons coming from other",40,pt,150.);
00179   Z_pt_postSelection = fs->make<TH1F>("Z_pt_postSelection","p_{T} distribution for muons coming from Z after iso selection",40,pt,150.);
00180   W_pt_postSelection = fs->make<TH1F>("W_pt_postSelection","p_{t} distribution for muons coming from W after iso selection",40,pt,150.);
00181   Other_pt_postSelection = fs->make<TH1F>("Other_pt_postSelection","p_{t} distribution for muons coming from other after iso selection",40,pt,150.);
00182 
00183 
00184   TkrPt = fs->make<TH1F>("TkrPt","IsoDeposit p distribution in the Tracker",100,0.,10.);
00185   EcalEt = fs->make<TH1F>("EcalEt","IsoDeposit E distribution in the Ecal",100,0.,5.);
00186   HcalEt = fs->make<TH1F>("HcalEt","IsoDeposit E distribution in the Hcal",100,0.,5.);
00187 
00188   EcalEtZ = fs->make<TH1F>("VetoEcalEt"," #Sigma E_{T} deposited in veto cone in the Ecal",100,0.,10.);
00189   HcalEtZ = fs->make<TH1F>("VetoHcalEt"," #Sigma E_{T} deposited in veto cone in the Hcal",100,0.,10.);
00190 }
00191 
00192 void ZMuMuIsolationAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
00193   Handle<CandidateView> dimuons;
00194   event.getByLabel(src_muons,dimuons);
00195 
00196   for(unsigned int i=0; i < dimuons->size(); ++ i ) {
00197     const Candidate & zmm = (* dimuons)[i];
00198     const Candidate * dau0 = zmm.daughter(0);
00199     const Candidate * dau1 = zmm.daughter(1);
00200     const pat::Muon & mu0 = dynamic_cast<const pat::Muon &>(*dau0->masterClone());
00201     const pat::GenericParticle & mu1 = dynamic_cast<const pat::GenericParticle &>(*dau1->masterClone());
00202 
00203     const pat::IsoDeposit * muTrackIso =mu0.isoDeposit(pat::TrackIso);
00204     const pat::IsoDeposit * tkTrackIso =mu1.isoDeposit(pat::TrackIso);
00205     const pat::IsoDeposit * muEcalIso =mu0.isoDeposit(pat::EcalIso);
00206     const pat::IsoDeposit * tkEcalIso =mu1.isoDeposit(pat::EcalIso);
00207     const pat::IsoDeposit * muHcalIso =mu0.isoDeposit(pat::HcalIso);
00208     const pat::IsoDeposit * tkHcalIso =mu1.isoDeposit(pat::HcalIso);
00209 
00210    
00211     if(mu0.pt() > pt && mu1.pt() > pt && abs(mu0.eta()) < eta && abs(mu1.eta()) < eta){
00212 
00213       Direction muDir = Direction(mu0.eta(),mu0.phi());
00214       Direction tkDir = Direction(mu1.eta(),mu1.phi());
00215 
00216       IsoDeposit::AbsVetos vetos_mu;
00217       vetos_mu.push_back(new ConeVeto( muDir, dRVeto ));
00218       vetos_mu.push_back(new ThresholdVeto( ptThreshold ));
00219 
00220       reco::IsoDeposit::AbsVetos vetos_tk;
00221       vetos_tk.push_back(new ConeVeto( tkDir, dRVeto ));
00222       vetos_tk.push_back(new ThresholdVeto( ptThreshold ));
00223       
00224       reco::IsoDeposit::AbsVetos vetos_mu_ecal;
00225       vetos_mu_ecal.push_back(new ConeVeto( muDir, 0. ));
00226       vetos_mu_ecal.push_back(new ThresholdVeto( etEcalThreshold ));
00227       
00228       reco::IsoDeposit::AbsVetos vetos_tk_ecal;
00229       vetos_tk_ecal.push_back(new ConeVeto( tkDir, 0. ));
00230       vetos_tk_ecal.push_back(new ThresholdVeto( etEcalThreshold ));
00231       
00232       reco::IsoDeposit::AbsVetos vetos_mu_hcal;
00233       vetos_mu_hcal.push_back(new ConeVeto( muDir, 0. ));
00234       vetos_mu_hcal.push_back(new ThresholdVeto( etHcalThreshold ));
00235       
00236       reco::IsoDeposit::AbsVetos vetos_tk_hcal;
00237       vetos_tk_hcal.push_back(new ConeVeto( tkDir, 0. ));
00238       vetos_tk_hcal.push_back(new ThresholdVeto( etHcalThreshold ));
00239       MuTag tag_mu = muTag(mu0);
00240       MuTag tag_track = muTag(mu1);
00241      
00242       double  Tk_isovalue = TMath::Max(muTrackIso->sumWithin(dRTrk,vetos_mu),tkTrackIso->sumWithin(dRTrk, vetos_tk));
00243       double  Ecal_isovalue = TMath::Max(muEcalIso->sumWithin(dREcal,vetos_mu_ecal),tkEcalIso->sumWithin(dREcal, vetos_tk_ecal));
00244       double  Hcal_isovalue = TMath::Max(muHcalIso->sumWithin(dRHcal,vetos_mu_hcal),tkHcalIso->sumWithin(dRHcal, vetos_tk_hcal));
00245       EcalEtZ->Fill(muEcalIso->candEnergy());
00246       EcalEtZ->Fill(tkEcalIso->candEnergy());
00247       HcalEtZ->Fill(muHcalIso->candEnergy());
00248       HcalEtZ->Fill(tkHcalIso->candEnergy());
00249 
00250       double iso_value0 = alpha*((0.5*(1+beta)* muEcalIso->sumWithin(dREcal,vetos_mu_ecal) ) + (0.5*(1-beta)*muHcalIso->sumWithin(dRHcal,vetos_mu_hcal) ) ) +(1-alpha)*muTrackIso->sumWithin(dRTrk,vetos_mu);
00251       double iso_value1 = alpha*((0.5*(1+beta)* tkEcalIso->sumWithin(dREcal,vetos_tk_ecal) ) + (0.5*(1-beta)*tkHcalIso->sumWithin(dRHcal,vetos_tk_hcal) ) ) +(1-alpha)*tkTrackIso->sumWithin(dRTrk,vetos_tk);
00252 
00253       double iso_value=TMath::Max(iso_value0,iso_value1);
00254 
00255       if(tag_mu==muFromZ && tag_track==muFromZ){
00256         h_IsoZ_tk->Fill(Tk_isovalue);
00257         h_IsoZ_ecal->Fill(Ecal_isovalue);
00258         h_IsoZ_hcal->Fill(Hcal_isovalue);
00259         IsoZ->Fill(iso_value);
00260 
00261         Z_eta->Fill(mu0.eta());
00262         Z_eta->Fill(mu1.eta());
00263         Z_pt->Fill(mu0.pt());
00264         Z_pt->Fill(mu1.pt());
00265 
00266         if(iso_value0<iso_cut)  {
00267           Z_pt_postSelection->Fill(mu0.pt());
00268           Z_eta_postSelection->Fill(mu0.eta());
00269         }
00270         if(iso_value1<iso_cut){
00271           Z_pt_postSelection->Fill(mu1.pt());
00272           Z_eta_postSelection->Fill(mu1.eta());
00273         }
00274 
00275         Deposits(muTrackIso,dRTrk,TkrPt);
00276         Deposits(muEcalIso,dREcal,EcalEt);
00277         Deposits(muHcalIso,dRHcal,HcalEt);
00278         Deposits(tkTrackIso,dRTrk,TkrPt);
00279         Deposits(tkEcalIso,dREcal,EcalEt);
00280         Deposits(tkHcalIso,dRHcal,HcalEt);
00281       }
00282       if(tag_mu==muFromW || tag_track==muFromW){ 
00283         h_IsoW_tk->Fill(Tk_isovalue);
00284         h_IsoW_ecal->Fill(Ecal_isovalue);
00285         h_IsoW_hcal->Fill(Hcal_isovalue);
00286         IsoW->Fill(iso_value);
00287 
00288         W_eta->Fill(mu0.eta());
00289         W_eta->Fill(mu1.eta());
00290         W_pt->Fill(mu0.pt());
00291         W_pt->Fill(mu1.pt());
00292 
00293         if(iso_value0<iso_cut) {
00294           W_pt_postSelection->Fill(mu0.pt());
00295           W_eta_postSelection->Fill(mu0.eta());
00296         }
00297         if(iso_value1<iso_cut) {
00298           W_pt_postSelection->Fill(mu1.pt());
00299           W_eta_postSelection->Fill(mu1.eta());
00300         }
00301 
00302         Deposits(muTrackIso,dRTrk,TkrPt);
00303         Deposits(muEcalIso,dREcal,EcalEt);
00304         Deposits(muHcalIso,dRHcal,HcalEt);
00305         Deposits(tkTrackIso,dRTrk,TkrPt);
00306         Deposits(tkEcalIso,dREcal,EcalEt);
00307         Deposits(tkHcalIso,dRHcal,HcalEt);
00308       }
00309       else{
00310         h_IsoOther_tk->Fill(Tk_isovalue);
00311         h_IsoOther_ecal->Fill(Ecal_isovalue);
00312         h_IsoOther_hcal->Fill(Hcal_isovalue);
00313         IsoOther->Fill(iso_value);
00314 
00315         Other_eta->Fill(mu0.eta());
00316         Other_eta->Fill(mu1.eta());
00317         Other_pt->Fill(mu0.pt());
00318         Other_pt->Fill(mu1.pt());
00319 
00320         if(iso_value0<iso_cut) {
00321           Other_pt_postSelection->Fill(mu0.pt());
00322           Other_eta_postSelection->Fill(mu0.eta());
00323         }
00324         if(iso_value1<iso_cut) {
00325           Other_pt_postSelection->Fill(mu1.pt());
00326           Other_eta_postSelection->Fill(mu1.eta());
00327         }
00328 
00329         Deposits(muTrackIso,dRTrk,TkrPt);
00330         Deposits(muEcalIso,dREcal,EcalEt);
00331         Deposits(muHcalIso,dRHcal,HcalEt);
00332         Deposits(tkTrackIso,dRTrk,TkrPt);
00333         Deposits(tkEcalIso,dREcal,EcalEt);
00334         Deposits(tkHcalIso,dRHcal,HcalEt);
00335       }
00336     }
00337   }    
00338 
00339   histo(h_IsoZ_tk,"#Sigma p_{T}","Events");
00340   histo(h_IsoW_tk,"#Sigma p_{T}","Events");
00341   histo(h_IsoOther_tk,"#Sigma p_{T}","#Events");
00342   histo(h_IsoZ_ecal,"#Sigma E_{t}","Events");
00343   histo(h_IsoW_ecal,"#Sigma E_{t}","Events");
00344   histo(h_IsoOther_ecal,"#Sigma E_{t}","Events");
00345   histo(h_IsoZ_hcal,"#Sigma E_{t}","Events");
00346   histo(h_IsoW_hcal,"#Sigma E_{t}","Events");
00347   histo(h_IsoOther_hcal,"#Sigma E_{t}","Events");
00348   histo(TkrPt,"p ","");
00349   histo(EcalEt,"E ","");
00350   histo(HcalEt,"E ","");
00351   histo(HcalEtZ,"E_{T}","");
00352   histo(EcalEtZ,"E_{T}","");
00353 }
00354   
00355 void ZMuMuIsolationAnalyzer::endJob() {
00356 }
00357 
00358 #include "FWCore/Framework/interface/MakerMacros.h"
00359 
00360 DEFINE_FWK_MODULE(ZMuMuIsolationAnalyzer);
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMuIsolationAnalyzer.cc
