New histograms for fast check of track reconstruction in the task to QA MC productions
authorfprino <fprino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Oct 2013 21:57:05 +0000 (21:57 +0000)
committerfprino <fprino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Oct 2013 21:57:05 +0000 (21:57 +0000)
PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx
PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h
PWGHF/vertexingHF/macros/PlotOutputMCCheck.C

index d5b4df8..af8b47a 100644 (file)
@@ -2,6 +2,7 @@
 #include "AliAnalysisManager.h"
 #include "AliAnalysisDataContainer.h"
 #include "AliESDEvent.h"
+#include "AliESDtrackCuts.h"
 #include "AliStack.h"
 #include "AliCentrality.h"
 #include "AliMCEventHandler.h"
@@ -17,6 +18,7 @@
 #include <TNtuple.h>
 #include <TH1F.h>
 #include <TH2F.h>
+#include <TH3F.h>
 #include <TChain.h>
 #include "AliESDInputHandlerRP.h"
 #include "AliAnalysisTaskCheckHFMCProd.h"
@@ -76,8 +78,17 @@ AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE
   fHistDSpecies(0),
   fHistBSpecies(0),
   fHistNcollHFtype(0),
+  fHistEtaPhiPtGenEle(0),
+  fHistEtaPhiPtGenPi(0),
+  fHistEtaPhiPtGenK(0),
+  fHistEtaPhiPtGenPro(0),
+  fHistEtaPhiPtRecEle(0),
+  fHistEtaPhiPtRecPi(0),
+  fHistEtaPhiPtRecK(0),
+  fHistEtaPhiPtRecPro(0),
   fSearchUpToQuark(kFALSE),
   fSystem(0),
+  fESDtrackCuts(0x0),
   fReadMC(kTRUE)
 {
   //
@@ -94,6 +105,8 @@ AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
     delete fOutput;
     fOutput = 0;
   }
+  delete fESDtrackCuts;
+
 }
    
 //___________________________________________________________________________
@@ -294,6 +307,39 @@ void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
 
   fHistNcollHFtype=new TH2F("hNcollHFtype","",5,-1.5,3.5,30,-0.5,29.5);
   fOutput->Add(fHistNcollHFtype);
+
+  Double_t binseta[11]={-1.0,-0.8,-0.6,-0.4,-0.2,0.,0.2,0.4,0.6,0.8,1.0};
+  const Int_t nBinsPhi=40;
+  Double_t binsphi[nBinsPhi+1];
+  for(Int_t ib=0; ib<=nBinsPhi; ib++) binsphi[ib]=ib*TMath::Pi()/20.;
+  const Int_t nBinsPt=24;  
+  Double_t binspt[nBinsPt+1]={0.,0.10,0.15,0.2,0.25,
+                             0.3,0.4,0.5,0.6,0.7,
+                             0.8,0.9,1.,1.25,1.5,
+                             1.75,2.,2.5,3.,4.,
+                             5.,7.5,10.,15.,20.};
+
+   fHistEtaPhiPtGenEle=new TH3F("hEtaPhiPtGenEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtGenEle);
+  fHistEtaPhiPtGenPi=new TH3F("hEtaPhiPtGenPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtGenPi);
+  fHistEtaPhiPtGenK=new TH3F("hEtaPhiPtGenK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtGenK);
+  fHistEtaPhiPtGenPro=new TH3F("hEtaPhiPtGenPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtGenPro);
+
+
+  fHistEtaPhiPtRecEle=new TH3F("hEtaPhiPtRecEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtRecEle);
+  fHistEtaPhiPtRecPi=new TH3F("hEtaPhiPtRecPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtRecPi);
+  fHistEtaPhiPtRecK=new TH3F("hEtaPhiPtRecK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtRecK);
+  fHistEtaPhiPtRecPro=new TH3F("hEtaPhiPtRecPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
+  fOutput->Add(fHistEtaPhiPtRecPro);
+
+
   PostData(1,fOutput);
 
 }
@@ -312,6 +358,18 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
 
   fHistoNEvents->Fill(0);
 
+  if(!fESDtrackCuts){
+    Int_t year=2011;
+    if(esd->GetRunNumber()<=139517) year=2010;
+    if(year==2010) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);
+    else fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE); 
+    fESDtrackCuts->SetMaxDCAToVertexXY(2.4);
+    fESDtrackCuts->SetMaxDCAToVertexZ(3.2);
+    fESDtrackCuts->SetDCAToVertex2D(kTRUE);
+    fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+                                           AliESDtrackCuts::kAny);
+  }
+
   Int_t nTracks=esd->GetNumberOfTracks();
   fHistoTracks->Fill(nTracks);
   Int_t nSelTracks=0;
@@ -374,6 +432,8 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
       Printf("ERROR: generated vertex not available");
       return;
     }
+    if(TMath::Abs(mcVert->GetZ())>10) return;
+
     //    const AliHeader* h=(AliHeader*)mcEvent->GetHeader();
     //    cout<<h<<endl;
     TString genname=mcEvent->GenEventHeader()->ClassName();
@@ -392,7 +452,6 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
        else if(title.Contains("cele")) typeHF=2;
       }
       nColl=lgen->GetEntries();
-      printf("Ncoll=%d typeHF=%d\n",nColl,typeHF);
       fHistNcollHFtype->Fill(typeHF,nColl);
     }
     Int_t nParticles=stack->GetNtrack();
@@ -409,6 +468,11 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
       if(stack->IsPhysicalPrimary(i)){
        Double_t eta=part->Eta();
        fHistoEtaPhysPrim->Fill(eta);
+       if(absPdg==11) fHistEtaPhiPtGenEle->Fill(eta,part->Phi(),part->Pt());
+       else if(absPdg==211) fHistEtaPhiPtGenPi->Fill(eta,part->Phi(),part->Pt());
+       else if(absPdg==321) fHistEtaPhiPtGenK->Fill(eta,part->Phi(),part->Pt());
+       else if(absPdg==2212) fHistEtaPhiPtGenPro->Fill(eta,part->Phi(),part->Pt());
+       
        if(TMath::Abs(eta)<0.5){
          dNchdy+=0.6666;   // 2/3 for the ratio charged/all
          nPhysPrim++;
@@ -550,6 +614,21 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
       else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);      
     }
 
+    for(Int_t i=0; i<nTracks; i++){
+      AliESDtrack* track=esd->GetTrack(i);
+      if(fESDtrackCuts->AcceptTrack(track)){
+       Int_t label=TMath::Abs(track->GetLabel());
+       if(stack->IsPhysicalPrimary(label)){
+         TParticle* part = (TParticle*)stack->Particle(label);
+         Int_t absPdg=TMath::Abs(part->GetPdgCode());
+         Double_t eta=part->Eta();
+         if(absPdg==11) fHistEtaPhiPtRecEle->Fill(eta,part->Phi(),part->Pt());
+         else if(absPdg==211) fHistEtaPhiPtRecPi->Fill(eta,part->Phi(),part->Pt());
+         else if(absPdg==321) fHistEtaPhiPtRecK->Fill(eta,part->Phi(),part->Pt());
+         else if(absPdg==2212) fHistEtaPhiPtRecPro->Fill(eta,part->Phi(),part->Pt());      
+       }
+      }
+    }
     fHistoNcharmed->Fill(dNchdy,nCharmed);
     fHistoNbVsNc->Fill(nc,nb);
     fHistoPhysPrim->Fill(nPhysPrim);
index eca40f5..66758b4 100644 (file)
@@ -19,11 +19,13 @@ class TList;
 class TNtuple;
 class TH1F;
 class TH2F;
+class TH3F;
 class TTree;
 class TString;
 class AliESDEvent;
 class AliESDfriend;
 class AliStack;
+class AliESDtrackCuts;
 
 #include "AliAnalysisTaskSE.h"
 
@@ -97,11 +99,20 @@ class AliAnalysisTaskCheckHFMCProd : public AliAnalysisTaskSE {
   TH1F* fHistDSpecies;          //! histo of D hadron species
   TH1F* fHistBSpecies;          //! histo of B hadron species
   TH2F* fHistNcollHFtype;      //! histo of B hadron species
+  TH3F* fHistEtaPhiPtGenEle;   //! histo of generated electrons
+  TH3F* fHistEtaPhiPtGenPi;   //! histo of generated pions
+  TH3F* fHistEtaPhiPtGenK;   //! histo of generated kaons
+  TH3F* fHistEtaPhiPtGenPro;   //! histo of generated protons
+  TH3F* fHistEtaPhiPtRecEle;   //! histo of generated electrons
+  TH3F* fHistEtaPhiPtRecPi;   //! histo of generated pions
+  TH3F* fHistEtaPhiPtRecK;   //! histo of generated kaons
+  TH3F* fHistEtaPhiPtRecPro;   //! histo of generated protons
   Bool_t fSearchUpToQuark; // c/b separation using quarks
   Int_t fSystem;         // 0=pp, 1=PbPb, 2=pPb
+  AliESDtrackCuts *fESDtrackCuts; // track selection
   Bool_t fReadMC;
 
-  ClassDef(AliAnalysisTaskCheckHFMCProd,5);
+  ClassDef(AliAnalysisTaskCheckHFMCProd,6);
 };
 
 
index 4bb600a..2fdc328 100644 (file)
 
 
 void PlotOutputMCCheck(){
+
   TFile *fil=new TFile("AnalysisResults.root");
   TDirectoryFile* df=(TDirectoryFile*)fil->Get("HFMCCheck");
   TList* l=(TList*)df->Get("clistHFMCCheck");
-  l->ls();
 
   TH1F* hNEvents=(TH1F*)l->FindObject("hNEvents");
   Int_t nAnalEv=hNEvents->GetBinContent(1);
   printf("Number of events= %d\n",nAnalEv);
 
-
   TCanvas* cv=new TCanvas("cv","Vertex");
   cv->Divide(3,3);
   cv->cd(1);
@@ -72,6 +71,243 @@ void PlotOutputMCCheck(){
   TH1F* hSelTracks=(TH1F*)l->FindObject("hSelTracks");
   hSelTracks->Draw();
 
+  // tracking efficiency
+  Double_t minEta=-0.8;
+  Double_t maxEta=0.8;
+  Double_t minPt=0.5;
+  Double_t maxPt=2.;
+  TH3F* hEtaPhiPtGenPi=(TH3F*)l->FindObject("hEtaPhiPtGenPi");
+  Int_t minEtaBin=hEtaPhiPtGenPi->GetXaxis()->FindBin(minEta+0.00001);
+  Int_t maxEtaBin=hEtaPhiPtGenPi->GetXaxis()->FindBin(maxEta-0.00001);
+  cout<<minEtaBin<<"    "<<maxEtaBin<<endl;
+  Int_t minPtBin=hEtaPhiPtGenPi->GetZaxis()->FindBin(minPt+0.00001);
+  Int_t maxPtBin=hEtaPhiPtGenPi->GetZaxis()->FindBin(maxPt-0.00001);
+  cout<<minPtBin<<"    "<<maxPtBin<<endl;
+  TH1D* hEtaGenPi=hEtaPhiPtGenPi->ProjectionX("hEtaGenPi",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiGenPi=hEtaPhiPtGenPi->ProjectionY("hPhiGenPi",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtGenPi=hEtaPhiPtGenPi->ProjectionZ("hPtGenPi",minEtaBin,maxEtaBin);
+  hEtaGenPi->Sumw2();
+  hPhiGenPi->Sumw2();
+  hPtGenPi->Sumw2();
+  hEtaGenPi->GetXaxis()->SetTitle("#eta");
+  hPhiGenPi->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtGenPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH3F* hEtaPhiPtRecPi=(TH3F*)l->FindObject("hEtaPhiPtRecPi");
+  TH1D* hEtaRecPi=hEtaPhiPtRecPi->ProjectionX("hEtaRecPi",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiRecPi=hEtaPhiPtRecPi->ProjectionY("hPhiRecPi",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtRecPi=hEtaPhiPtRecPi->ProjectionZ("hPtRecPi",minEtaBin,maxEtaBin);
+  hEtaRecPi->Sumw2();
+  hPhiRecPi->Sumw2();
+  hPtRecPi->Sumw2();
+  hEtaRecPi->GetXaxis()->SetTitle("#eta");
+  hPhiRecPi->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtRecPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH1D* hEtaEffPi=(TH1D*)hEtaRecPi->Clone("hEtaEffPi");
+  hEtaEffPi->Divide(hEtaRecPi,hEtaGenPi,1.,1.,"B");
+  TH1D* hPhiEffPi=(TH1D*)hPhiRecPi->Clone("hPhiEffPi");
+  hPhiEffPi->Divide(hPhiRecPi,hPhiGenPi,1.,1.,"B");
+  TH1D* hPtEffPi=(TH1D*)hPtRecPi->Clone("hPtEffPi");
+  hPtEffPi->Divide(hPtRecPi,hPtGenPi,1.,1.,"B");
+
+  TH3F* hEtaPhiPtGenK=(TH3F*)l->FindObject("hEtaPhiPtGenK");
+  TH1D* hEtaGenK=hEtaPhiPtGenK->ProjectionX("hEtaGenK",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiGenK=hEtaPhiPtGenK->ProjectionY("hPhiGenK",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtGenK=hEtaPhiPtGenK->ProjectionZ("hPtGenK",minEtaBin,maxEtaBin);
+  hEtaGenK->Sumw2();
+  hPhiGenK->Sumw2();
+  hPtGenK->Sumw2();
+  hEtaGenK->GetXaxis()->SetTitle("#eta");
+  hPhiGenK->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtGenK->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH3F* hEtaPhiPtRecK=(TH3F*)l->FindObject("hEtaPhiPtRecK");
+  TH1D* hEtaRecK=hEtaPhiPtRecK->ProjectionX("hEtaRecK",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiRecK=hEtaPhiPtRecK->ProjectionY("hPhiRecK",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtRecK=hEtaPhiPtRecK->ProjectionZ("hPtRecK",minEtaBin,maxEtaBin);
+  hEtaRecK->Sumw2();
+  hPhiRecK->Sumw2();
+  hPtRecK->Sumw2();
+  hEtaRecK->GetXaxis()->SetTitle("#eta");
+  hPhiRecK->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtRecK->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH1D* hEtaEffK=(TH1D*)hEtaRecK->Clone("hEtaEffK");
+  hEtaEffK->Divide(hEtaRecK,hEtaGenK,1.,1.,"B");
+  TH1D* hPhiEffK=(TH1D*)hPhiRecK->Clone("hPhiEffK");
+  hPhiEffK->Divide(hPhiRecK,hPhiGenK,1.,1.,"B");
+  TH1D* hPtEffK=(TH1D*)hPtRecK->Clone("hPtEffK");
+  hPtEffK->Divide(hPtRecK,hPtGenK,1.,1.,"B");
+
+  TH3F* hEtaPhiPtGenPro=(TH3F*)l->FindObject("hEtaPhiPtGenPro");
+  TH1D* hEtaGenPro=hEtaPhiPtGenPro->ProjectionX("hEtaGenPro",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiGenPro=hEtaPhiPtGenPro->ProjectionY("hPhiGenPro",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtGenPro=hEtaPhiPtGenPro->ProjectionZ("hPtGenPro",minEtaBin,maxEtaBin);
+  hEtaGenPro->Sumw2();
+  hPhiGenPro->Sumw2();
+  hPtGenPro->Sumw2();
+  hEtaGenPro->GetXaxis()->SetTitle("#eta");
+  hPhiGenPro->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtGenPro->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH3F* hEtaPhiPtRecPro=(TH3F*)l->FindObject("hEtaPhiPtRecPro");
+  TH1D* hEtaRecPro=hEtaPhiPtRecPro->ProjectionX("hEtaRecPro",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiRecPro=hEtaPhiPtRecPro->ProjectionY("hPhiRecPro",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtRecPro=hEtaPhiPtRecPro->ProjectionZ("hPtRecPro",minEtaBin,maxEtaBin);
+  hEtaRecPro->Sumw2();
+  hPhiRecPro->Sumw2();
+  hPtRecPro->Sumw2();
+  hEtaRecPro->GetXaxis()->SetTitle("#eta");
+  hPhiRecPro->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtRecPro->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH1D* hEtaEffPro=(TH1D*)hEtaRecPro->Clone("hEtaEffPro");
+  hEtaEffPro->Divide(hEtaRecPro,hEtaGenPro,1.,1.,"B");
+  TH1D* hPhiEffPro=(TH1D*)hPhiRecPro->Clone("hPhiEffPro");
+  hPhiEffPro->Divide(hPhiRecPro,hPhiGenPro,1.,1.,"B");
+  TH1D* hPtEffPro=(TH1D*)hPtRecPro->Clone("hPtEffPro");
+  hPtEffPro->Divide(hPtRecPro,hPtGenPro,1.,1.,"B");
+
+  TH3F* hEtaPhiPtGenEle=(TH3F*)l->FindObject("hEtaPhiPtGenEle");
+  TH1D* hEtaGenEle=hEtaPhiPtGenEle->ProjectionX("hEtaGenEle",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiGenEle=hEtaPhiPtGenEle->ProjectionY("hPhiGenEle",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtGenEle=hEtaPhiPtGenEle->ProjectionZ("hPtGenEle",minEtaBin,maxEtaBin);
+  hEtaGenEle->Sumw2();
+  hPhiGenEle->Sumw2();
+  hPtGenEle->Sumw2();
+  hEtaGenEle->GetXaxis()->SetTitle("#eta");
+  hPhiGenEle->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtGenEle->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH3F* hEtaPhiPtRecEle=(TH3F*)l->FindObject("hEtaPhiPtRecEle");
+  TH1D* hEtaRecEle=hEtaPhiPtRecEle->ProjectionX("hEtaRecEle",0,-1,minPtBin,maxPtBin);
+  TH1D* hPhiRecEle=hEtaPhiPtRecEle->ProjectionY("hPhiRecEle",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
+  TH1D* hPtRecEle=hEtaPhiPtRecEle->ProjectionZ("hPtRecEle",minEtaBin,maxEtaBin);
+  hEtaRecEle->Sumw2();
+  hPhiRecEle->Sumw2();
+  hPtRecEle->Sumw2();
+  hEtaRecEle->GetXaxis()->SetTitle("#eta");
+  hPhiRecEle->GetXaxis()->SetTitle("#varphi (rad)");
+  hPtRecEle->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  TH1D* hEtaEffEle=(TH1D*)hEtaRecEle->Clone("hEtaEffEle");
+  hEtaEffEle->Divide(hEtaRecEle,hEtaGenEle,1.,1.,"B");
+  TH1D* hPhiEffEle=(TH1D*)hPhiRecEle->Clone("hPhiEffEle");
+  hPhiEffEle->Divide(hPhiRecEle,hPhiGenEle,1.,1.,"B");
+  TH1D* hPtEffEle=(TH1D*)hPtRecEle->Clone("hPtEffEle");
+  hPtEffEle->Divide(hPtRecEle,hPtGenEle,1.,1.,"B");
+
+  for(Int_t iBin=1; iBin<=hPtGenPi->GetNbinsX(); iBin++){
+    Double_t w=hPtGenPi->GetBinWidth(iBin);
+    Double_t c=hPtGenPi->GetBinContent(iBin);
+    hPtGenPi->SetBinContent(iBin,c/w);
+    c=hPtGenK->GetBinContent(iBin);
+    hPtGenK->SetBinContent(iBin,c/w);
+    c=hPtGenPro->GetBinContent(iBin);
+    hPtGenPro->SetBinContent(iBin,c/w);
+    c=hPtGenEle->GetBinContent(iBin);
+    hPtGenEle->SetBinContent(iBin,c/w);
+
+  }
+
+  hEtaEffPi->SetStats(0);
+  hPtEffPi->SetStats(0);
+  hPhiEffPi->SetStats(0);
+  hEtaEffPi->SetMinimum(0.);
+  hPtEffPi->SetMinimum(0.);
+  hPhiEffPi->SetMinimum(0.);
+  hEtaEffPi->SetMaximum(1.05);
+  hPtEffPi->SetMaximum(1.05);
+  hPhiEffPi->SetMaximum(1.05);
+
+  hPtGenPi->SetMarkerStyle(27);
+  hPtEffPi->SetMarkerStyle(27);
+  hEtaEffPi->SetMarkerStyle(27);
+  hPhiEffPi->SetMarkerStyle(27);
+  hPtGenK->SetMarkerStyle(22);
+  hPtEffK->SetMarkerStyle(22);
+  hEtaEffK->SetMarkerStyle(22);
+  hPhiEffK->SetMarkerStyle(22);
+  hPtGenK->SetMarkerColor(2);
+  hPtEffK->SetMarkerColor(2);
+  hEtaEffK->SetMarkerColor(2);
+  hPhiEffK->SetMarkerColor(2);
+  hPtGenPro->SetMarkerStyle(25);
+  hPtEffPro->SetMarkerStyle(25);
+  hEtaEffPro->SetMarkerStyle(25);
+  hPhiEffPro->SetMarkerStyle(25);
+  hPtGenPro->SetMarkerColor(4);
+  hPtEffPro->SetMarkerColor(4);
+  hEtaEffPro->SetMarkerColor(4);
+  hPhiEffPro->SetMarkerColor(4);
+  hPtGenEle->SetMarkerStyle(20);
+  hPtEffEle->SetMarkerStyle(20);
+  hEtaEffEle->SetMarkerStyle(20);
+  hPhiEffEle->SetMarkerStyle(20);
+  hPtGenEle->SetMarkerColor(kGreen+1);
+  hPtEffEle->SetMarkerColor(kGreen+1);
+  hEtaEffEle->SetMarkerColor(kGreen+1);
+  hPhiEffEle->SetMarkerColor(kGreen+1);
+
+  TCanvas* ctref=new TCanvas("ctref","Track eff",1200,800);
+  ctref->Divide(2,2);
+  ctref->cd(1);
+  gPad->SetLogy();
+  hPtGenPi->Draw();
+  ctref->Update();
+  TPaveStats* st1=(TPaveStats*)hPtGenPi->GetListOfFunctions()->FindObject("stats");
+  st1->SetY2NDC(0.91);
+  st1->SetY1NDC(0.76);
+  st1->SetTextColor(1);
+  hPtGenK->Draw("sames");
+  ctref->Update();
+  TPaveStats* st2=(TPaveStats*)hPtGenK->GetListOfFunctions()->FindObject("stats");
+  st2->SetY2NDC(0.75);
+  st2->SetY1NDC(0.60);
+  st2->SetTextColor(2);
+  hPtGenPro->Draw("sames");
+  ctref->Update();
+  TPaveStats* st3=(TPaveStats*)hPtGenPro->GetListOfFunctions()->FindObject("stats");
+  st3->SetY2NDC(0.59);
+  st3->SetY1NDC(0.44);
+  st3->SetTextColor(4);
+  hPtGenEle->Draw("sames");
+  ctref->Update();
+  TPaveStats* st4=(TPaveStats*)hPtGenEle->GetListOfFunctions()->FindObject("stats");
+  st4->SetY2NDC(0.43);
+  st4->SetY1NDC(0.28);
+  st4->SetTextColor(kGreen+1);
+  gPad->Modified();
+  TLegend* leg=new TLegend(0.5,0.5,0.7,0.8);
+  leg->SetFillColor(0);
+  leg->SetBorderSize(0);
+  leg->AddEntry(hPtGenPi,"Pions","P")->SetTextColor(hPtGenPi->GetMarkerColor());
+  leg->AddEntry(hPtGenK,"Kaons","P")->SetTextColor(hPtGenK->GetMarkerColor());
+  leg->AddEntry(hPtGenPro,"Protons","P")->SetTextColor(hPtGenPro->GetMarkerColor());
+  leg->AddEntry(hPtGenEle,"Electrons","P")->SetTextColor(hPtGenEle->GetMarkerColor());
+  leg->Draw();
+
+  gPad->Update();
+  ctref->cd(2);
+  hPtEffPi->Draw();
+  hPtEffPi->GetYaxis()->SetTitle("Efficiency");
+  hPtEffK->Draw("same");
+  hPtEffEle->Draw("same");
+  hPtEffPro->Draw("same");
+  TLatex* t1=new TLatex(0.62,0.18,Form("%.2f < #eta < %.2f",minEta,maxEta));
+  t1->SetNDC();
+  t1->Draw();
+  ctref->cd(3);
+  hPhiEffPi->Draw();
+  hPhiEffPi->GetYaxis()->SetTitle("Efficiency");
+  hPhiEffK->Draw("same");
+  hPhiEffEle->Draw("same");
+  hPhiEffPro->Draw("same");
+  TLatex* t2=new TLatex(0.18,0.18,Form("%.1f < p_{T} < %.1f GeV/c",minPt,maxPt));
+  t2->SetNDC();
+  t2->Draw();
+  t1->Draw();
+  ctref->cd(4);
+  hEtaEffPi->Draw();
+  hEtaEffPi->GetYaxis()->SetTitle("Efficiency");
+  hEtaEffK->Draw("same");
+  hEtaEffEle->Draw("same");
+  hEtaEffPro->Draw("same");
+  t2->Draw();
+
   TH1F* hncharmed=(TH1F*)l->FindObject("hncharmed");
   TCanvas* cn=new TCanvas("cn","ncharm");
   hncharmed->Draw("box");