PWGLF/totEt/macros/hadEt/CorrEfficiency.C

   1 //Christine Nattrass, University of Tennessee at Knoxville
   2 //This macro is for calculating the single track efficiency for the hadronic transverse energy measurement
   3 //Uses the output of AliAnalysisTaskHadEt
   4 //This is not actually what gets used in the correction class AliAnalysisHadEtCorrections - that is done in the macro GetCorrections.C - but this is useful for making plots and playing around with different options
   5
   6 TFile *file;// = new TFile(filename);
   7
   8 //this function calculates the efficiency propagating errors properly
   9 TH1D* bayneseffdiv(TH1D* numerator, TH1D* denominator,Char_t* name)
  10 {
  11     if(!numerator){
  12       cerr<<"Error:  numerator does not exist!"<<endl;
  13       return NULL;
  14     }
  15     if(!denominator){
  16       cerr<<"Error:  denominator does not exist!"<<endl;
  17       return NULL;
  18     }
  19     TH1D* result = (TH1D*)numerator->Clone(name);
  20     Int_t nbins = numerator->GetNbinsX();
  21     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
  22       Double_t numeratorVal = numerator->GetBinContent(ibin);
  23       Double_t denominatorVal = denominator->GetBinContent(ibin);
  24       // Check if the errors are right or the thing is scaled
  25       Double_t numeratorValErr = numerator->GetBinError(ibin);
  26       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
  27         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
  28         numeratorVal /= rescale;
  29       }
  30       Double_t denominatorValErr = denominator->GetBinError(ibin);
  31       if (!(denominatorValErr==0. || denominatorVal==0. )) {
  32         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
  33         denominatorVal /= rescale;
  34       }
  35       Double_t quotient = 0.;
  36       if (denominatorVal!=0.) {
  37         quotient = numeratorVal/denominatorVal;
  38       }
  39       Double_t quotientErr=0;
  40       quotientErr = TMath::Sqrt(
  41                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
  42                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
  43       result->SetBinContent(ibin,quotient);
  44       result->SetBinError(ibin,quotientErr);
  45       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
  46     }
  47     return result;
  48 }
  49
  50
  51 //This is a somewhat messy function that gets the efficiency for different particles
  52 TH1D *GetHisto(float cut = 0.12, char *name, int mycase, bool eta, int color, int marker,bool TPC, bool ITS, bool PbPb, int cb, int cblast, char *filename){
  53   //TFile *file = new TFile("Et.ESD.new.sim.merged.root");
  54   //TFile *file = new TFile("rootFiles/LHC11a4_bis/Et.ESD.new.sim.LHC11a4_bis.root");
  55   file->cd();
  56   //TFile *file = new TFile("rootFiles/LHC10d7/Et.ESD.new.sim.LHC10d7.root");
  57   TList *list = file->FindObject("out2");
  58   char *myname = "ITS";
  59   if(TPC&&!ITS){ delete myname; myname = "TPC";}
  60   if(TPC&&ITS) {delete myname; myname = "TPCITS";}
  61   char *mynamealt = "TPCITS";
  62   //cout<<"Using tracks from "<<myname<<" for efficiency"<<endl;
  63   TH2F *numeratorParent;
  64   switch(mycase){
  65   case 0:
  66     if(cblast != -1){//add more centrality bins
  67       for(int i=cb;i<=cblast;i++){
  68         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)))->Clone("RecoHadron");
  69         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)));}
  70         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiMinus",i)));
  71         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KMinus",i)));
  72         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KPlus",i)));
  73         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"Proton",i)));
  74         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"AntiProton",i)));
  75         if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
  76           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"PiPlus",i)));
  77           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"PiMinus",i)));
  78           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"KMinus",i)));
  79           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"KPlus",i)));
  80           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"Proton",i)));
  81           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"AntiProton",i)));
  82         }
  83       }
  84     }
  85     else{
  86       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoHadron");
  87       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
  88       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
  89       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")));
  90       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")));
  91       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
  92       if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
  93         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"PiPlus")));
  94         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"PiMinus")));
  95         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"KMinus")));
  96         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"KPlus")));
  97         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"Proton")));
  98         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"AntiProton")));
  99       }
 100     }
 101     break;
 102   case 1://pion
 103     if(cblast != -1){//add more centrality bins
 104       for(int i=cb;i<=cblast;i++){
 105         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)))->Clone("RecoPion");
 106         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiPlus",i)));}
 107         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"PiMinus",i)));
 108         if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
 109           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"PiMinus",i)));
 110           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"PiPlus",i)));
 111         }
 112       }
 113     }
 114     else{
 115       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiPlus")))->Clone("RecoPion");
 116       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"PiMinus")));
 117       if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
 118         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"PiPlus")));
 119         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"PiMinus")));
 120       }
 121     }
 122     break;
 123   case 2://kaon
 124     if(cblast != -1){//add more centrality bins
 125       for(int i=cb;i<=cblast;i++){
 126         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KPlus",i)))->Clone("RecoKaon");
 127         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KPlus",i)));}
 128         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"KMinus",i)));
 129         if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
 130           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"KPlus",i)));
 131           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"KMinus",i)));
 132         }
 133       }
 134     }
 135     else{
 136       // cout<<"I am kaoning here !"<<endl;
 137       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KPlus")))->Clone("RecoKaon");
 138       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"KMinus")));
 139       if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
 140         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"KPlus")));
 141         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"KMinus")));
 142       }
 143       //cout<<"Numerator "<<numeratorParent->GetEntries()<<endl;
 144     }
 145     break;
 146   case 3://proton
 147     if(cblast != -1){//add more centrality bins
 148       for(int i=cb;i<=cblast;i++){
 149         if(i==cb) numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"Proton",i)))->Clone("RecoProton");
 150         else{numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"Proton",i)));}
 151         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",myname,"AntiProton",i)));
 152         if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
 153           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"Proton",i)));
 154           numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%sCB%i",mynamealt,"AntiProton",i)));
 155         }
 156       }
 157     }
 158     else{
 159       //cout<<"I am protoning here !"<<endl;
 160       numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"Proton")))->Clone("RecoProton");
 161       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",myname,"AntiProton")));
 162       if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
 163         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"Proton")));
 164         numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s",mynamealt,"AntiProton")));
 165       }
 166       //cout<<"Numerator "<<numeratorParent->GetEntries()<<endl;
 167     }
 168     break;
 169   case 4://electron
 170     numeratorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s%s",myname,"EPlus",cbname)))->Clone("RecoElectron");
 171     numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s%s",myname,"EMinus",cbname)));
 172     if(!TPC&&ITS){//ITS Standalone tracks - from leftover hits, don't want to double count
 173       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s%s",mynamealt,"EPlus",cbname)));
 174       numeratorParent->Add((TH2F*) out2->FindObject(Form("EtNReconstructed%s%s%s",mynamealt,"EMinus",cbname)));
 175     }
 176     break;
 177   }
 178   TH2F *denominatorParent;
 179   switch(mycase){
 180   case 0:
 181     if(cblast != -1){//add more centrality bins
 182       denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedChargedHadronCB%i",cb)))->Clone("RecoHadron");
 183       for(int i=cb+1;i<=cblast;i++){
 184         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedChargedHadronCB%i",i)));
 185       }
 186     }
 187     else{
 188       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedChargedHadron"))->Clone("RecoHadron");
 189     }
 190     break;
 191   case 1://pion
 192     if(cblast != -1){//add more centrality bins
 193       denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedPiPlusCB%i",cb)))->Clone("RecoPion");
 194       denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedPiMinusCB%i",cb)));
 195       for(int i=cb+1;i<=cblast;i++){
 196         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedPiPlusCB%i",i)));
 197         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedPiMinusCB%i",i)));
 198       }
 199     }
 200     else{
 201       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedPiPlus"))->Clone("RecoPion");
 202       denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedPiMinus"));
 203     }
 204     break;
 205   case 2://kaon
 206     if(cblast != -1){//add more centrality bins
 207       denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedKPlusCB%i",cb)))->Clone("RecoKaon");
 208       denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedKMinusCB%i",cb)));
 209       for(int i=cb+1;i<=cblast;i++){
 210         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedKPlusCB%i",i)));
 211         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedKMinusCB%i",i)));
 212       }
 213     }
 214     else{
 215       //cout<<"I am here kaoning"<<endl;
 216       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedKPlus"))->Clone("RecoKaon");
 217       denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedKMinus"));
 218       //cout<<"Denominator "<<denominatorParent->GetEntries()<<endl;
 219     }
 220     break;
 221   case 3://proton
 222     if(cblast != -1){//add more centrality bins
 223       for(int i=cb;i<=cblast;i++){
 224         if(cb==i)denominatorParent = (TH2F*)((TH2F*) out2->FindObject(Form("EtNSimulatedProtonCB%i",i)))->Clone("RecoProton");
 225         else{denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedProtonCB%i",i)));}
 226         denominatorParent->Add((TH2F*) out2->FindObject(Form("EtNSimulatedAntiProtonCB%i",i)));
 227       }
 228     }
 229     else{
 230       //cout<<"I am here protoning"<<endl;
 231       denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedProton"))->Clone("RecoProton");
 232       denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedAntiProton"));
 233       //cout<<"Denominator "<<denominatorParent->GetEntries()<<endl;
 234     }
 235     break;
 236   case 4://electron
 237     denominatorParent = (TH2F*)((TH2F*) out2->FindObject("EtNSimulatedEPlus"))->Clone("RecoElectron");
 238     denominatorParent->Add((TH2F*) out2->FindObject("EtNSimulatedEMinus"));
 239     break;
 240   }
 241   //cout<<"Numerator "<<numeratorParent->GetEntries()<<endl;
 242   // cout<<"Denominator "<<denominatorParent->GetEntries()<<endl;
 243   numeratorParent->Sumw2();
 244   denominatorParent->Sumw2();
 245   TH1D *denominator;
 246   TH1D *numerator;
 247   if(eta){
 248     int lowbin = numeratorParent->GetYaxis()->FindBin(-cut+.001);//make sure we don't accv0entally get the wrong bin
 249     int highbin = numeratorParent->GetYaxis()->FindBin(cut-.001);
 250     //cout<<"Projecting from "<<numeratorParent->GetYaxis()->GetBinLowEdge(lowbin)<<" to "<<numeratorParent->GetYaxis()->GetBinLowEdge(highbin+1)<<endl;
 251     denominator = denominatorParent->ProjectionX(Form("garbage%s",name),lowbin,highbin);
 252     numerator = numeratorParent->ProjectionX(name,lowbin,highbin);
 253   }
 254   else{
 255     int lowbin = denominatorParent->GetXaxis()->FindBin(cut);//make sure we don't accidentally get the wrong bin
 256     int highbin = denominatorParent->GetXaxis()->GetNbins();
 257     // cout<<"Here Projecting from "<<denominatorParent->GetXaxis()->GetBinLowEdge(lowbin)<<" to "<<denominatorParent->GetXaxis()->GetBinLowEdge(highbin+1)<<endl;
 258     numerator = numeratorParent->ProjectionY(name,lowbin,highbin);
 259     denominator = denominatorParent->ProjectionY(Form("denominator%s",name),lowbin,highbin);
 260   }
 261   delete numeratorParent;
 262   delete denominatorParent;
 263   //numerator->Divide(denominator);
 264   TH1D *result = bayneseffdiv((TH1D*) numerator,(TH1D*)denominator,name);
 265   //result->Rebin(2);
 266   //result->Scale(0.5);
 267   result->SetYTitle("Efficiency");
 268   result->GetYaxis()->SetTitleOffset(0.8);
 269   result->GetXaxis()->SetTitleOffset(0.8);
 270   result->GetYaxis()->SetLabelSize(0.05);
 271   result->GetXaxis()->SetLabelSize(0.05);
 272   result->GetYaxis()->SetTitleSize(0.05);
 273   result->GetXaxis()->SetTitleSize(0.05);
 274   result->SetMarkerColor(color);
 275   result->SetLineColor(color);
 276   result->SetMarkerStyle(marker);
 277   //result->Draw("e");
 278   //file->Close();
 279   //delete file;
 280   delete numerator;
 281   delete denominator;
 282   result->SetName(name);
 283   return result;
 284
 285 }
 286
 287
 288 //this is a method that makes pretty plots
 289 void CorrEfficiency(char *prodname = "LHC11a4_bis HIJING 2.76 TeV Pb+Pb",char *shortprodname= "LHC11a4_bis", char *filename = "rootFiles/LHC11a4_bis/Et.ESD.new.sim.LHC11a4_bis.root", bool TPC = false,bool ITS = true, bool eta = true, bool PbPb = true, int cb = 0, int cblast = -1){
 290
 291   file = new TFile(filename);
 292   gStyle->SetOptTitle(0);
 293   gStyle->SetOptStat(0);
 294   gStyle->SetOptFit(0);
 295   TCanvas *c = new TCanvas("c","c",600,400);
 296   c->SetTopMargin(0.02);
 297   c->SetRightMargin(0.02);
 298   c->SetBorderSize(0);
 299   c->SetFillColor(0);
 300   c->SetFillColor(0);
 301   c->SetBorderMode(0);
 302   c->SetFrameFillColor(0);
 303   c->SetFrameBorderMode(0);
 304   //c->SetLogx();
 305
 306   int colortotal = 1;
 307   int colorpi = 2;
 308   int colork = 3;
 309   int colorp = 4;
 310   int phosmarker = 20;
 311   int emcalmarker = 24;
 312   float ptcut1 = 0.05;
 313   float ptcut2 = 0.1;
 314   float phoscut = 0.12;
 315   float emcalcut = 0.7;
 316   if(!eta){
 317     phoscut = 0.1;
 318     emcalcut = 0.15;
 319   }
 320   TH1D *PHOStotal = GetHisto(phoscut,"PHOStotal",0,eta,colortotal,phosmarker,TPC,ITS,PbPb,cb,cblast,filename);
 321   TH1D *PHOSpi = GetHisto(phoscut,"PHOSpi",1,eta,colorpi,phosmarker,TPC,ITS,PbPb,cb,cblast,filename);
 322   TH1D *PHOSp = GetHisto(phoscut,"PHOSp",3,eta,colork,phosmarker,TPC,ITS,PbPb,cb,cblast,filename);
 323   TH1D *PHOSk = GetHisto(phoscut,"PHOSk",2,eta,colorp,phosmarker,TPC,ITS,PbPb,cb,cblast,filename);
 324   if(eta) PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));
 325 //if(ITS&&!TPC){PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.05),PHOStotal->GetXaxis()->FindBin(1.0));}
 326 //else{PHOStotal->GetXaxis()->SetRange(PHOStotal->GetXaxis()->FindBin(0.0),PHOStotal->GetXaxis()->FindBin(3.0));}
 327   PHOStotal->SetMinimum(0.0);
 328   PHOStotal->SetMaximum(1.0);
 329   //parameters[centbin][0]*exp(-pow(parameters[centbin][1]/pt,parameters[centbin][2]))
 330   TF1 *func = new TF1("func","[0]*exp(-pow([1]/x,[2]))");
 331   func->SetParameter(0,.9);
 332   func->SetParameter(1,.05);
 333   func->SetParLimits(1,1e-3,1);
 334   func->SetParameter(2,.1);
 335   //PHOStotal->Fit(func);
 336   PHOStotal->Draw();
 337   PHOSpi->Draw("same");
 338   PHOSp->Draw("same");
 339   PHOSk->Draw("same");
 340   //cout<<"Hadrons"<<endl;
 341   TH1D *EMCALtotal = GetHisto(emcalcut,"EMCALtotal",0,eta,colortotal,emcalmarker,TPC,ITS,PbPb,cb,cblast,filename);
 342   //cout<<endl<<endl<<"=================================PIONS================================="<<endl;
 343   TH1D *EMCALpi = GetHisto(emcalcut,"EMCALpi",1,eta,colorpi,emcalmarker,TPC,ITS,PbPb,cb,cblast,filename);
 344   //cout<<endl<<endl<<"=================================PROTONS================================="<<endl;
 345   TH1D *EMCALp = GetHisto(emcalcut,"EMCALp",3,eta,colork,emcalmarker,TPC,ITS,PbPb,cb,cblast,filename);
 346   //cout<<endl<<endl<<"=================================KAONS================================="<<endl;
 347   TH1D *EMCALk = GetHisto(emcalcut,"EMCALk",2,eta,colorp,emcalmarker,TPC,ITS,PbPb,cb,cblast,filename);
 348   EMCALtotal->Draw("same");
 349   EMCALpi->Draw("same");
 350   EMCALp->Draw("same");
 351   EMCALk->Draw("same");
 352
 353
 354   TLegend *leg = new  TLegend(0.22651,0.247312,0.370805,0.438172);
 355   leg->AddEntry(PHOStotal,"#pi,K,p");
 356   leg->AddEntry(PHOSpi,"#pi^{#pm}");
 357   leg->AddEntry(PHOSk,"K^{#pm}");
 358   leg->AddEntry(PHOSp,"p,#bar{p}");
 359   leg->SetFillStyle(0);
 360   leg->SetFillColor(0);
 361   leg->SetBorderSize(0);
 362   leg->SetTextSize(0.06);
 363  leg->Draw();
 364
 365   TLine *line = new TLine(0.150,0.0,0.150,1.0);
 366   if(eta)line->Draw();
 367   line->SetLineWidth(3.0);
 368   //line->SetLineColor(TColor::kYellow);
 369   line->SetLineStyle(2);
 370   TLatex *tex = new TLatex(0.398954,0.0513196,prodname);
 371   tex->SetTextSize(0.0537634);
 372   tex->Draw();
 373   TLatex *tex3 = new TLatex(1.16186,0.28348,"Closed symbols |#eta|<0.12 (PHOS)");
 374   tex3->SetTextSize(0.0537634);
 375   tex3->Draw();
 376   TLatex *tex4 = new TLatex(1.16186,0.213221,"Open symbols |#eta|<0.70 (EMCal)");
 377   tex4->SetTextSize(0.0537634);
 378   tex4->Draw();
 379   TLatex *tex2 = new TLatex(0.164016,0.860826,"TPC cut-off 150 MeV/c");
 380   tex2->SetTextSize(0.0537634);
 381   if(eta) tex2->Draw();
 382
 383
 384   TLine *line2 = new TLine(0.10,0.0,0.10,1.0);
 385   line2->SetLineWidth(3.0);
 386   TLatex *tex5 = new TLatex(0.10817,0.924976,"ITS cut-off 100 MeV/c");
 387   tex5->SetTextSize(0.0537634);
 388   line2->SetLineStyle(2);
 389   tex5->SetTextColor(4);
 390   line2->SetLineColor(4);
 391   if(!TPC && eta){
 392     line2->Draw();
 393     tex5->Draw();
 394   }
 395   if(!TPC){
 396     c->SaveAs(Form("pics/%s/CorrEfficiency.eps",shortprodname));
 397     c->SaveAs(Form("pics/%s/CorrEfficiency.png",shortprodname));
 398   }
 399   else{
 400     c->SaveAs(Form("pics/%s/CorrEfficiencyTPCITS.eps",shortprodname));
 401     c->SaveAs(Form("pics/%s/CorrEfficiencyTPCITS.png",shortprodname));
 402   }
 403
 404   if(PbPb){//make one more plot
 405     //pions - no real centrality dependence
 406     //three centrality bins for efficiency 0-25%, 25-50%, 50-90%
 407     //use same for unidentified hadrons
 408     //kaons & protons - centrality dependence is more significant but I don't think I can do better on the binning
 409     //int pid = 2;//h=0,pi=1,p=3,k=2
 410     int maxid = 3;
 411     if(!TPC) maxid=0;
 412     for(int pid=0;pid<=maxid;pid++){
 413       TCanvas *c2 = new TCanvas("c2","c2",600,400);
 414       c2->SetTopMargin(0.02);
 415       c2->SetRightMargin(0.02);
 416       c2->SetBorderSize(0);
 417       c2->SetFillColor(0);
 418       c2->SetFillColor(0);
 419       c2->SetBorderMode(0);
 420       c2->SetFrameFillColor(0);
 421       c2->SetFrameBorderMode(0);
 422       //cout<<endl<<endl;
 423
 424       TH1D *cb0 = GetHisto(phoscut,"cb0",pid,eta,1,20,TPC,ITS,PbPb,0,4,filename);
 425       TH1D *cb4 = GetHisto(phoscut,"cb5",pid,eta,4,20,TPC,ITS,PbPb,5,9,filename);
 426       //TH1D *cb9 = GetHisto(phoscut,"cb9",pid,eta,TColor::kGreen+4,20,TPC,ITS,PbPb,10,14);
 427       TH1D *cb14 = GetHisto(phoscut,"cb14",pid,eta,2,20,TPC,ITS,PbPb,10,18,filename);
 428       cb0->GetXaxis()->SetRange(cb0->GetXaxis()->FindBin(0.05),cb0->GetXaxis()->FindBin(5.0));
 429       cb0->SetMaximum(1.0);
 430       cb0->Draw();
 431       cb4->Draw("same");
 432       // cb9->Draw("same");
 433       cb14->Draw("same");
 434       TLegend *leg = new TLegend(0.124161,0.747312,0.318792,0.959677);//(left,bottom,right,top)
 435       leg->SetTextSize(0.0537634);
 436       leg->SetBorderSize(0);
 437       leg->SetLineColor(0);
 438       leg->SetLineStyle(0);
 439       leg->SetLineWidth(0);
 440       leg->SetFillColor(0);
 441       leg->SetFillStyle(0);
 442       leg->AddEntry(cb0,"0-25%");
 443       leg->AddEntry(cb4,"25-50%");
 444       //leg->AddEntry(cb9,"45-50%");
 445       leg->AddEntry(cb14,"50-90%");
 446       leg->Draw();
 447       TString particle;
 448       TString unidentified = "h";
 449       TString proton = "p";
 450       TString pion = "#pi";
 451       TString kaon = "K";
 452       switch(pid){
 453       case 0:
 454         particle = unidentified;
 455         break;
 456       case 1:
 457         particle = pion;
 458         break;
 459       case 2:
 460         particle = kaon;
 461         break;
 462       case 3:
 463         particle = proton;
 464         break;
 465       }
 466       TLatex *texpid = new TLatex(0.25237,0.589541,particle.Data());
 467       texpid->SetTextSize(0.0806452);
 468       texpid->Draw();
 469       switch(pid){
 470       case 0:
 471         if(!TPC){c2->SaveAs(Form("pics/%s/CorrEfficiencyITSUnidentified.eps",shortprodname,pid));}
 472         else{c2->SaveAs(Form("pics/%s/CorrEfficiencyTPCITSUnidentified.eps",shortprodname,pid));}
 473         break;
 474       case 1:
 475         if(!TPC){c2->SaveAs(Form("pics/%s/CorrEfficiencyITSPion.eps",shortprodname,pid));}
 476         else{c2->SaveAs(Form("pics/%s/CorrEfficiencyTPCITSPion.eps",shortprodname,pid));}
 477         break;
 478       case 3:
 479         if(!TPC){c2->SaveAs(Form("pics/%s/CorrEfficiencyITSProton.eps",shortprodname,pid));}
 480         else{c2->SaveAs(Form("pics/%s/CorrEfficiencyTPCITSProton.eps",shortprodname,pid));}
 481         break;
 482       case 2:
 483         if(!TPC){c2->SaveAs(Form("pics/%s/CorrEfficiencyITSKaon.eps",shortprodname,pid));}
 484         else{c2->SaveAs(Form("pics/%s/CorrEfficiencyTPCITSKaon.eps",shortprodname,pid));}
 485         break;
 486       }
 487     }
 488   }
 489 }