PWGHF/vertexingHF/charmFlow/DmesonsFlowAnalysis.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TFile.h>
   4 #include <TString.h>
   5 #include <TH2F.h>
   6 #include <TH1F.h>
   7 #include <TF1.h>
   8 #include <TGraph.h>
   9 #include <TGraphErrors.h>
  10 #include <TMultiGraph.h>
  11 #include <TDirectoryFile.h>
  12 #include <TList.h>
  13 #include <TCanvas.h>
  14 #include <TLegend.h>
  15 #include <TLegendEntry.h>
  16 #include <TPaveText.h>
  17 #include <TStyle.h>
  18 #include <TASImage.h>
  19 #include <TPad.h>
  20 #include <TROOT.h>
  21 #include <TDatabasePDG.h>
  22 #include <TParameter.h>
  23 #include <TLatex.h>
  24 #include <TGraphAsymmErrors.h>
  25
  26 #include <AliHFMassFitter.h>
  27 #include "AliRDHFCutsD0toKpi.h"
  28 #include "AliVertexingHFUtils.h"
  29 #include "AliRDHFCutsDplustoKpipi.h"
  30 #include "AliRDHFCutsDStartoKpipi.h"
  31 #include <AliVertexingHFUtils.h>
  32
  33 #endif
  34
  35 //methods for the analysis of AliAnalysisTaskSEHFv2 output
  36 //Authors: Chiara Bianchin cbianchi@pd.infn.it
  37 //         Giacomo Ortona  ortona@to.infn.it
  38 //         Francesco Prino prino@to.infn.it
  39
  40 //global variables to be set
  41 Bool_t gnopng=kTRUE; //don't save in png format (only root and eps)
  42 const Int_t nptbinsnew=3;
  43 Float_t ptbinsnew[nptbinsnew+1]={3,5,8,12.};
  44 Int_t fittype=0;
  45 Int_t rebin[nptbinsnew]={4,4,4};
  46 Double_t nsigma=3;
  47 const Int_t nphibins=4;
  48 Float_t phibinslim[nphibins+1]={0,TMath::Pi()/4,TMath::Pi()/2,3*TMath::Pi()/4,TMath::Pi()};
  49 Int_t minPtBin[nptbinsnew]={-1,-1,-1};
  50 Int_t maxPtBin[nptbinsnew]={-1,-1,-1};
  51 Double_t mass;
  52 //methods
  53 //Bool_t ReadFile(TList* &list,TH1F* &hstat,AliRDHFCuts* &cutobj,TString listname,TString partname,TString path="./",TString filename="AnalysisResults.root");
  54 TList *LoadMassHistos(TList *inputlist,Int_t minCent,Int_t maxCent,Bool_t inoutanis);
  55 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value);
  56 void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs,Bool_t inoutanis);
  57 void DmesonsFlowAnalysis(Bool_t inoutanis=kTRUE,Int_t minC=30,Int_t maxC=50,TString partname="Dplus",Int_t evPlane=2);
  58 void DrawEventPlane(Int_t mincentr=0,Int_t maxcentr=0,TString filename="AnalysisResults.root",TString dirname="PWG3_D2H_HFv2",TString listname="coutputv2",Int_t evPlane=2);
  59 Double_t DrawEventPlane(TList *list,Int_t mincentr=0,Int_t maxcentr=0,Int_t evPlane=2);
  60 void DrawEventPlaneAllCentralities(TList *list,Int_t startCentr=20,Int_t endCentr=80,Int_t evPlane=2);
  61 //Aggiungere <pt> method
  62
  63 //methods implementation
  64 //________________________________________________________________________________
  65 TList *LoadMassHistos(TList *inputlist,Int_t minCent,Int_t maxCent,Bool_t inoutanis){
  66   // printf("Start load histos\n");
  67   //  const Int_t nptbins=cutobj->GetNPtBins();
  68   TList *outlist = new TList();
  69   outlist->SetName("azimuthalhistoslist");
  70
  71   Int_t nphi=nphibins;
  72   if(inoutanis)nphi=2;
  73
  74   //Create 2D histogram in final pt bins
  75   for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
  76     for(Int_t iphi=0;iphi<nphi;iphi++){
  77       TH1F *hMass=0x0;//=new TH1F();
  78       for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin];iPtBin++){
  79         for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){
  80           TString hisname=Form("hMphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+5);
  81           TH2F* htmp=(TH2F*)inputlist->FindObject(hisname.Data());
  82           Int_t startX=htmp->FindBin(phibinslim[iphi]);
  83           Int_t endX=htmp->FindBin(phibinslim[iphi+1]);
  84           TH1F *h1tmp;
  85           if(inoutanis){
  86             if(iphi==0){
  87               h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi0",iPtBin),htmp->FindBin(0),htmp->FindBin(TMath::Pi()/4.));
  88               h1tmp->Add((TH1F*)htmp->ProjectionY(Form("hMass%d",iPtBin),htmp->FindBin(3.*TMath::Pi()/4.),htmp->FindBin(TMath::Pi())));
  89             }else{
  90               h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi1",iPtBin),htmp->FindBin(TMath::Pi()/4.),htmp->FindBin(3.*TMath::Pi()/4.));
  91             }
  92           }else h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi%d",iPtBin,iphi),startX,endX);
  93           if(hMass==0)hMass=(TH1F*)h1tmp->Clone();
  94           else hMass->Add((TH1F*)h1tmp->Clone());
  95         }
  96       }
  97       hMass->SetTitle(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
  98       hMass->SetName(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
  99       outlist->Add(hMass->Clone());
 100       //      hMass->DrawClone();
 101       delete hMass;
 102       hMass=0x0;
 103     }
 104   }
 105   return outlist;
 106 }
 107 //______________________________________________________________
 108 Bool_t DefinePtBins(AliRDHFCuts *cutobj){
 109   Int_t nPtBinsCuts=cutobj->GetNPtBins();
 110   Float_t *ptlimsCuts=cutobj->GetPtBinLimits();
 111   //  for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f-%f\n",iPt,ptlimsCuts[iPt],ptlimsCuts[iPt+1]);
 112   for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
 113     for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
 114       if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[iFinalPtBin])<0.0001){
 115         minPtBin[iFinalPtBin]=iPtCuts;
 116         if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts-1;
 117       }
 118     }
 119     if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[nptbinsnew])<0.0001) maxPtBin[nptbinsnew-1]=iPtCuts-1;
 120   }
 121   if(TMath::Abs(ptbinsnew[nptbinsnew]-ptlimsCuts[nPtBinsCuts])<0.0001) maxPtBin[nptbinsnew-1]=nPtBinsCuts-1;
 122   for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
 123     printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
 124     if(minPtBin[iFinalPtBin]<0 || maxPtBin[iFinalPtBin]<0) return kFALSE;
 125   }
 126
 127   return kTRUE;
 128 }
 129 //______________________________________________________________
 130 Int_t GetPadNumber(Int_t ix,Int_t iy){
 131   return (iy)*nptbinsnew+ix+1;
 132 }
 133 //________________________________________________________________________________
 134 void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs,Bool_t inoutanis){
 135
 136   Int_t nphi=nphibins;
 137   if(inoutanis)nphi=2;
 138
 139   //Canvases for drawing histograms
 140   TCanvas *cDeltaPhi = new TCanvas("cinvmassdeltaphi","Invariant mass distributions",1200,700);
 141   TCanvas *cDeltaPhifs = new TCanvas("cinvmassdeltaphifs","Invariant mass distributions - fit with fixed sigma",1200,700);
 142   TCanvas *cPhiInteg = new TCanvas("cinvmass","Invariant mass distributions - #phi integrated",1200,350);
 143   cDeltaPhi->Divide(nptbinsnew,nphi);
 144   cDeltaPhifs->Divide(nptbinsnew,nphi);
 145   cPhiInteg->Divide(nptbinsnew,1);
 146
 147   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
 148     TH1F *histtofitfullsigma=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi0",ipt))->Clone();
 149     for(Int_t iphi=0;iphi<nphi;iphi++){
 150       Int_t ipad=GetPadNumber(ipt,iphi);
 151       Double_t signal=0,esignal=0;
 152       TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
 153       if(iphi>0)histtofitfullsigma->Add((TH1F*)histtofit->Clone());
 154       if(!histtofit){
 155         gSignal[ipt]->SetPoint(iphi,iphi,signal);
 156         gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
 157         return;
 158       }
 159       histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
 160       AliHFMassFitter fitter(histtofit,histtofit->GetBinLowEdge(2),histtofit->GetBinLowEdge(histtofit->GetNbinsX()-2),rebin[ipt]);
 161       fitter.SetInitialGaussianMean(mass);
 162       //      fitter.SetInitialGaussianSigma(0.012);
 163       Bool_t ok=fitter.MassFitter(kFALSE);
 164       if(ok){
 165         fitter.DrawHere(cDeltaPhi->cd(ipad),3,1);
 166         fitter.Signal(3,signal,esignal);
 167       }
 168       gSignal[ipt]->SetPoint(iphi,iphi,signal);
 169       gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
 170     }
 171     //fit for fixed sigma
 172     histtofitfullsigma->SetTitle(Form("%.1f<p_{t}<%.1f",ptbinsnew[ipt],ptbinsnew[ipt+1]));
 173     AliHFMassFitter fitter(histtofitfullsigma,histtofitfullsigma->GetBinLowEdge(2),histtofitfullsigma->GetBinLowEdge(histtofitfullsigma->GetNbinsX()-2),rebin[ipt]);
 174     fitter.SetInitialGaussianMean(mass);
 175     //      fitter.SetInitialGaussianSigma(0.012);
 176     Bool_t ok=fitter.MassFitter(kFALSE);
 177     if(ok){
 178       fitter.DrawHere(cPhiInteg->cd(ipt+1),3,1);
 179     }
 180     Double_t sigma=fitter.GetSigma();
 181     for(Int_t iphi=0;iphi<nphi;iphi++){
 182       Int_t ipad=GetPadNumber(ipt,iphi);
 183       TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
 184       histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
 185       AliHFMassFitter fitter2(histtofit,histtofit->GetBinLowEdge(2),histtofit->GetBinLowEdge(histtofit->GetNbinsX()-2),rebin[ipt]);
 186       fitter2.SetInitialGaussianMean(mass);
 187       fitter2.SetFixGaussianSigma(sigma);
 188       Bool_t ok2=fitter2.MassFitter(kFALSE);
 189       Double_t signal=0,esignal=0;
 190       if(ok2){
 191         fitter2.DrawHere(cDeltaPhifs->cd(ipad),3,1);
 192         fitter2.Signal(3,signal,esignal);
 193       }
 194       gSignalfs[ipt]->SetPoint(iphi,iphi,signal);
 195       gSignalfs[ipt]->SetPointError(iphi,0,0,esignal,esignal);
 196
 197     }
 198   }//end loop on pt bin
 199
 200   cDeltaPhi->SaveAs("InvMassDeltaPhi.eps");
 201   cDeltaPhifs->SaveAs("InvMassDeltaPhi_fs.eps");
 202   cPhiInteg->SaveAs("InvMassfullphi.eps");
 203
 204   if(!gnopng){
 205     cDeltaPhi->SaveAs("InvMassDeltaPhi.png");
 206     cDeltaPhifs->SaveAs("InvMassDeltaPhi_fs.png");
 207     cPhiInteg->SaveAs("InvMassfullphi.png");
 208   }
 209   //cDeltaPhifs->DrawClone();
 210   cDeltaPhifs->Close();
 211
 212 }
 213 //______________________________________________________________
 214 void DmesonsFlowAnalysis(Bool_t inoutanis,Int_t minC,Int_t maxC,TString partname,Int_t evPlane){
 215   TString filename="AnalysisResults.root";
 216   TString dirname="PWG3_D2H_HFv2";
 217   TString listname="coutputv2";
 218
 219   AliRDHFCuts *cutsobj=0x0;
 220   //Load input data from AliAnalysisTaskSEHFv2
 221   TFile *f = TFile::Open(filename.Data());
 222   if(!f){
 223     printf("file %s not found, please check file name\n",filename.Data());return;
 224   }
 225   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
 226   if(!dir){
 227     printf("Directory %s not found, please check dir name\n",dirname.Data());return;
 228   }
 229   if(partname.Contains("D0")) {
 230     cutsobj=((AliRDHFCutsD0toKpi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 231     mass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 232   }
 233   if(partname.Contains("Dplus")){
 234     cutsobj=((AliRDHFCutsDplustoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 235     mass=TDatabasePDG::Instance()->GetParticle(411)->Mass();
 236   }
 237   if(partname.Contains("Dstar")) {
 238     cutsobj=((AliRDHFCutsDStartoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
 239     mass=(TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass());
 240   }
 241
 242   TList *list =(TList*)dir->Get(listname.Data());
 243   if(!list){
 244     printf("list %s not found in file, please check list name\n",listname.Data());return;
 245   }
 246   if(!cutsobj){
 247     printf("cut object not found in file, please check keylist number\n");return;
 248   }
 249   //Define new pt bins
 250   if(!DefinePtBins(cutsobj)){
 251     printf("cut not define pt bins\n");return;
 252   }
 253   //Load mass histograms corresponding to the required centrality, pt range and phi binning
 254   printf("Load mass histos \n");
 255   TList *histlist=LoadMassHistos(list,minC,maxC,inoutanis);
 256   TString aniss="";
 257   if(inoutanis)aniss+="anis";
 258   histlist->SaveAs(Form("v2Histograms_%d_%d_%s.root",minC,maxC,aniss.Data()),"RECREATE");
 259
 260   Int_t nphi=nphibins;
 261   if(inoutanis)nphi=2;
 262   Int_t nSubRes=1;
 263   if(evPlane<2)nSubRes=3;//3 sub-events method for VZERO EP
 264
 265   //EP resolution
 266   TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",minC,minC+5));
 267   TH1F* hevplresos2=0;
 268   TH1F* hevplresos3=0;
 269   if(evPlane<2){
 270     hevplresos2=(TH1F*)list->FindObject(Form("hEvPlaneReso2centr%d_%d",minC,minC+5));
 271     hevplresos3=(TH1F*)list->FindObject(Form("hEvPlaneReso3centr%d_%d",minC,minC+5));
 272   }
 273   Double_t resol=0;
 274   for(Int_t icent=minC+5;icent<maxC;icent=icent+5){
 275     hevplresos->Add((TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",icent,icent+5)));
 276     if(evPlane<2){
 277       hevplresos2->Add((TH1F*)list->FindObject(Form("hEvPlaneReso2centr%d_%d",icent,icent+5)));
 278       hevplresos3->Add((TH1F*)list->FindObject(Form("hEvPlaneReso3centr%d_%d",icent,icent+5)));
 279     }
 280   }
 281   if(evPlane>=2){
 282     resol=AliVertexingHFUtils::GetFullEvResol(hevplresos);
 283   }else{
 284     Double_t resolSub[3]={hevplresos->GetMean(),hevplresos2->GetMean(),hevplresos3->GetMean()};
 285     resol=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
 286   }
 287
 288   printf("average pt for pt bin \n");
 289   //average pt for pt bin
 290   AliVertexingHFUtils *utils=new AliVertexingHFUtils();
 291   TH2F* hmasspt=(TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",minC,minC+5));
 292   for(Int_t icent=minC+5;icent<maxC;icent=icent+5)hmasspt->Add((TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",icent,icent+5)));
 293   Float_t averagePt[nptbinsnew];
 294   Float_t errorPt[nptbinsnew];
 295   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
 296     TH1F *histtofit = (TH1F*)hmasspt->ProjectionY("_py",hmasspt->FindBin(ptbinsnew[ipt]),hmasspt->FindBin(ptbinsnew[ipt+1]));
 297     AliHFMassFitter fitter(histtofit,histtofit->GetBinLowEdge(2),histtofit->GetBinLowEdge(histtofit->GetNbinsX()-2),1);
 298     fitter.MassFitter(kFALSE);
 299     Float_t massFromFit=fitter.GetMean();
 300     Float_t sigmaFromFit=fitter.GetSigma();
 301     TF1* funcB2=fitter.GetBackgroundRecalcFunc();
 302     utils->AveragePt(averagePt[ipt],errorPt[ipt],ptbinsnew[ipt],ptbinsnew[ipt+1],hmasspt,massFromFit,sigmaFromFit,funcB2);
 303   }
 304
 305   printf("Fill TGraphs for signal \n");
 306   //Fill TGraphs for signal
 307   TGraphAsymmErrors *gSignal[nptbinsnew];
 308   TGraphAsymmErrors *gSignalfs[nptbinsnew];
 309   for(Int_t i=0;i<nptbinsnew;i++){
 310     gSignal[i]=new TGraphAsymmErrors(nphi);
 311     gSignal[i]->SetName(Form("gasigpt%d",i));
 312     gSignal[i]->SetTitle(Form("Signal %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
 313     gSignal[i]->SetMarkerStyle(25);
 314     gSignalfs[i]=new TGraphAsymmErrors(nphi);
 315     gSignalfs[i]->SetName(Form("gasigfspt%d",i));
 316     gSignalfs[i]->SetTitle(Form("Signal (fixed sigma) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
 317     gSignalfs[i]->SetMarkerStyle(21);
 318   }
 319   FillSignalGraph(histlist,gSignal,gSignalfs,inoutanis);
 320
 321   printf("Compute v2 \n");
 322   //compute v2
 323   TCanvas *cv2 =new TCanvas("v2","v2");
 324   TCanvas *cv2fs =new TCanvas("v2_fs","v2 Fit with fixed sigma");
 325   TGraphAsymmErrors *gv2=new TGraphAsymmErrors(nptbinsnew);
 326   TGraphAsymmErrors *gv2fs=new TGraphAsymmErrors(nptbinsnew);
 327   gv2->SetName("gav2");
 328   gv2->SetTitle("v_{2}(p_{t})");
 329   gv2fs->SetName("gav2fs");
 330   gv2fs->SetTitle("v_{2}(p_{t}) (fit with fixed sigma)");
 331
 332   //Prepare output file
 333   TFile *fout=new TFile(Form("v2Output_%d_%d_%s.root",minC,maxC,aniss.Data()),"RECREATE");
 334
 335   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
 336     fout->cd();
 337     gSignal[ipt]->Write();
 338     gSignalfs[ipt]->Write();
 339
 340     if(inoutanis){
 341       //v2 from in-out anisotropy
 342       Double_t *y,*yfs;
 343       y=gSignal[ipt]->GetY();
 344       yfs=gSignalfs[ipt]->GetY();
 345       Double_t nIn=y[0];
 346       Double_t nOut=y[1];
 347       Double_t enIn=gSignal[ipt]->GetErrorY(0);
 348       Double_t enOut=gSignal[ipt]->GetErrorY(1);
 349       Double_t anis=(nIn-nOut)/(nIn+nOut);
 350       Double_t eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
 351       Double_t v2=anis*TMath::Pi()/4./resol;
 352       Double_t ev2=eAnis*TMath::Pi()/4./resol;
 353       gv2->SetPoint(ipt,averagePt[ipt],v2);
 354       gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 355
 356       nIn=yfs[0];
 357       nOut=yfs[1];
 358       enIn=gSignal[ipt]->GetErrorY(0);
 359       enOut=gSignal[ipt]->GetErrorY(1);
 360       anis=(nIn-nOut)/(nIn+nOut);
 361       eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
 362       v2=anis*TMath::Pi()/4./resol;
 363       ev2=eAnis*TMath::Pi()/4./resol;
 364       gv2fs->SetPoint(ipt,averagePt[ipt],v2);
 365       gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 366     }else{
 367       TF1 *flowFunc = new TF1("flow","[0]*(1.+2.*[1]*TMath::Cos(2.*x))");
 368       //v2 from fit to Deltaphi distribution
 369       gSignal[ipt]->Fit(flowFunc);
 370       Double_t v2 = flowFunc->GetParameter(1)/resol;
 371       Double_t ev2=flowFunc->GetParError(1)/resol;
 372       gv2->SetPoint(ipt,averagePt[ipt],v2);
 373       gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 374
 375       gSignalfs[ipt]->Fit(flowFunc);
 376       v2 = flowFunc->GetParameter(1)/resol;
 377       ev2=flowFunc->GetParError(1)/resol;
 378       gv2fs->SetPoint(ipt,averagePt[ipt],v2);
 379       gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
 380     }
 381   }
 382   gv2->SetMarkerStyle(20);
 383   gv2fs->SetMarkerStyle(20);
 384   cv2->cd();
 385   gv2->Draw("AP");
 386   gv2->GetYaxis()->SetTitle("v_{2}");
 387   gv2->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 388   cv2fs->cd();
 389   gv2fs->Draw("AP");
 390   gv2fs->GetYaxis()->SetTitle("v_{2}");
 391   gv2fs->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 392
 393   fout->cd();
 394   gv2->Write();
 395   gv2fs->Write();
 396   printf("Event plane resolution %f\n",resol);
 397   printf("Average pt\n");
 398   for(Int_t ipt=0;ipt<nptbinsnew;ipt++) printf("%f +- %f\n",averagePt[ipt],errorPt[ipt]);
 399 }
 400 //_______________________________________________________________________________________________________________________________
 401 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value){
 402   for (Int_t i=0;i<nbins;i++){
 403     if(value>=array[i] && value<array[i+1]){
 404       return i;
 405     }
 406   }
 407   cout<<value<< " out of range "<<array[0]<<", "<<array[nbins]<<endl;
 408   return -1;
 409 }
 410
 411 //_______________________________________________________________________________________________________________________________
 412 void DrawEventPlane(Int_t mincentr,Int_t maxcentr,TString filename,TString dirname,TString listname,Int_t evPlane){
 413   TFile *f = TFile::Open(filename.Data());
 414   if(!f){
 415     printf("file %s not found, please check file name\n",filename.Data());return;
 416   }
 417   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
 418   if(!dir){
 419     printf("Directory %s not found, please check dir name\n",dirname.Data());return;
 420   }
 421   TList *list =(TList*)dir->Get(listname.Data());
 422   if(!list){
 423     printf("list %s not found in file, please check list name\n",listname.Data());return;
 424   }
 425   DrawEventPlane(list,mincentr,maxcentr,evPlane);
 426 }
 427 //_______________________________________________________________________________________________________________________________
 428 Double_t DrawEventPlane(TList *list,Int_t mincentr,Int_t maxcentr,Int_t evPlane){
 429
 430   //draw the histograms correlated to the event plane, returns event plane resolution
 431
 432   gStyle->SetCanvasColor(0);
 433   gStyle->SetTitleFillColor(0);
 434   gStyle->SetStatColor(0);
 435   //gStyle->SetPalette(1);
 436   gStyle->SetOptFit(1);
 437
 438   Double_t resolFull=0;
 439   Int_t nSubRes=1;
 440   if(evPlane<2)nSubRes=3;//3 sub-events method for VZERO EP
 441   TString namereso[3]={"Reso","Reso2","Reso3"};
 442   TString suffixcentr=Form("centr%d_%d",mincentr,maxcentr);
 443   TH2F* hevpls=(TH2F*)list->FindObject(Form("hEvPlanecentr%d_%d",mincentr,mincentr+5));
 444   hevpls->SetName(Form("hEvPlane%s",suffixcentr.Data()));
 445   hevpls->SetTitle(Form("Event Plane angle %s",suffixcentr.Data()));
 446     //    TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));QUI!!!
 447   TH1F* hevplresos[3];
 448   for(Int_t ires=0;ires<nSubRes;ires++){
 449     hevplresos[ires]=(TH1F*)list->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),mincentr,mincentr+5));
 450     hevplresos[ires]->SetName(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
 451     hevplresos[ires]->SetTitle(Form("Event Plane Resolution %s%s",namereso[ires].Data(),suffixcentr.Data()));
 452   }
 453
 454   for(Int_t icentr=mincentr+5;icentr<maxcentr;icentr=icentr+5){
 455     TH2F* h=(TH2F*)list->FindObject(Form("hEvPlanecentr%d_%d",icentr,icentr+5));
 456     if(h)hevpls->Add(h);
 457     else cout<<"skipping ev plane "<<icentr<<"_"<<icentr+5<<endl;
 458     for(Int_t ires=0;ires<nSubRes;ires++){
 459       TH1F *hr=(TH1F*)list->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),icentr,icentr+5));
 460       if(hr)hevplresos[ires]->Add(hr);
 461       else cout<<"skipping ev pl reso "<<icentr<<"_"<<icentr+5<<endl;
 462     }
 463   }
 464
 465   TCanvas* cvtotevpl=new TCanvas(Form("cvtotevpl%s",suffixcentr.Data()),Form("Ev plane %s",suffixcentr.Data()),1200,400);
 466   cvtotevpl->Divide(3,1);
 467   cvtotevpl->cd(1);
 468   hevpls->Draw("COLZ");
 469   TH1F* htpc = (TH1F*)hevpls->ProjectionX();
 470   cvtotevpl->cd(2);
 471   htpc->Draw();
 472   htpc->Fit("pol0");
 473   TH1F* hv0 = (TH1F*)hevpls->ProjectionY();
 474   cvtotevpl->cd(3);
 475   hv0->Draw();
 476   hv0->Fit("pol0");
 477
 478   TCanvas* cvtotevplreso=new TCanvas(Form("cvtotevplreso%s",suffixcentr.Data()),Form("Ev plane Resolution %s",suffixcentr.Data()));
 479   cvtotevplreso->cd();
 480   hevplresos[0]->SetTitle("TPC cos2(#Psi_{A}-#Psi_{B})");
 481   hevplresos[0]->Draw();
 482   gStyle->SetOptStat(0);
 483   if(nSubRes>1){
 484     hevplresos[1]->SetLineColor(2);
 485     hevplresos[2]->SetLineColor(3);
 486     hevplresos[0]->SetTitle("cos2(#Psi_{V0A}-#Psi_{V0B})");
 487     hevplresos[1]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0A})");
 488     hevplresos[2]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0B})");
 489     hevplresos[1]->Draw("SAME");
 490     hevplresos[2]->Draw("SAME");
 491   }
 492   TLegend *leg = cvtotevplreso->BuildLegend();
 493   leg->SetLineColor(0);
 494   leg->SetFillColor(0);
 495
 496   if(nSubRes<=1) resolFull=AliVertexingHFUtils::GetFullEvResol(hevplresos[0]);
 497   else{
 498     Double_t resolSub[3];
 499     for(Int_t ires=0;ires<nSubRes;ires++)resolSub[ires]=hevplresos[ires]->GetMean();
 500     resolFull=TMath::Sqrt(resolSub[2]*resolSub[0]/resolSub[1]);
 501   }
 502
 503   TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
 504   pvreso->SetBorderSize(0);
 505   pvreso->SetFillStyle(0);
 506   pvreso->AddText(Form("Resolution on full event = %.4f\n",resolFull));
 507   pvreso->Draw();
 508
 509   TFile* fout=new TFile(Form("EvPlanecentr%d-%d.root",mincentr,maxcentr),"recreate");
 510   fout->cd();
 511   hevpls->Write();
 512   for(Int_t ires=0;ires<nSubRes;ires++)hevplresos[ires]->Write();
 513
 514   return resolFull;
 515 }
 516
 517 //_______________________________________________________________________________________________________________________________
 518 void DrawEventPlaneAllCentralities(TList *list,Int_t startCentr,Int_t endCentr,Int_t evPlane){
 519   TGraph* gresovscentr=new TGraph(0);
 520   gresovscentr->SetName("gresovscentr");
 521   gresovscentr->SetTitle(Form("Resolution vs Centrality;centrality (%s);Resolution","%"));
 522
 523   Int_t iPoint=0;
 524   Int_t step=5;//5% centrality step
 525   //  Int_t nCC=(endCentr-startCentr)/step;
 526
 527   for(Int_t i=startCentr;i<endCentr;i=i+step){
 528     Double_t resolFull=DrawEventPlane(list,i,i+step,evPlane);
 529     gresovscentr->SetPoint(iPoint,i+(Float_t)step/2.,resolFull);
 530     iPoint++;
 531   }
 532
 533   TCanvas* cresovscentr=new TCanvas("cresovscentr", "Resolution vs Centrality");
 534   cresovscentr->cd();
 535   gresovscentr->SetMarkerStyle(20);
 536   gresovscentr->Draw("AP");
 537   if(!gnopng) gresovscentr->SaveAs(Form("%s.png",gresovscentr->GetName()));
 538   gresovscentr->SaveAs(Form("%s.eps",gresovscentr->GetName()));
 539   TFile* fout=new TFile(Form("%s.root",gresovscentr->GetName()),"recreate");
 540   fout->cd();
 541   gresovscentr->Write();
 542 }
 543 //_______________________________________________________________________________________________________________________________
 544 void DrawPaveStat(TVirtualPad* c,Int_t nev,Int_t mincentr,Int_t maxcentr,TString string,TString string2){
 545
 546   TPaveText* txtoff=new TPaveText(0.2,0.3,0.6,0.5,"NDC");
 547   txtoff->SetBorderSize(0);
 548   txtoff->SetFillStyle(0);
 549   //txtoff->AddText(Form("%.0f#times10^{6} events",nev/1e6));
 550   txtoff->AddText("Pb-Pb, #sqrt{s_{NN}} = 2.76TeV");
 551   printf("Centrality range %d%s, nev=%d, corrected=%f \n",maxcentr-mincentr,"%",nev,(Float_t)nev/(Float_t)(maxcentr-mincentr)/100.);
 552   txtoff->AddText(Form("%.0f#times10^{6} evts in %d-%d%s centr cl",nev/1e6*(Float_t)(maxcentr-mincentr)/100.,mincentr,maxcentr,"%"));
 553   txtoff->AddText(string.Data());
 554   if(string2!="") {
 555     txtoff->AddText(string2.Data());
 556     txtoff->SetY1(0.3);
 557   }
 558   txtoff->SetTextColor(kAzure-5);
 559   c->cd();
 560   txtoff->Draw();
 561 }
 562 //_______________________________________________________________________________________________________________________________
 563 void ALICEPerformance2(TVirtualPad *c,TString today,Double_t* where,TString type){
 564
 565   //date must be in the form: 04/05/2010
 566   if(today=="today"){
 567     TDatime startt;
 568     int date=startt.GetDate();
 569     int y=date/10000;
 570     int m=(date%10000)/100;
 571     int d=date%100;
 572
 573
 574     today="";
 575     today+=d;
 576     if(m<10)
 577       today.Append("/0");
 578     else today.Append("/");
 579     today+=m;
 580     today.Append("/");
 581     today+=y;
 582
 583   }
 584   cout<<"here"<<endl;
 585   TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",where[0],where[1],where[2],where[3]);
 586   myPadLogo->SetFillColor(0);
 587   myPadLogo->SetBorderMode(0);
 588   myPadLogo->SetBorderSize(2);
 589   myPadLogo->SetFrameBorderMode(0);
 590   myPadLogo->SetLeftMargin(0.0);
 591   myPadLogo->SetTopMargin(0.0);
 592   myPadLogo->SetBottomMargin(0.0);
 593   myPadLogo->SetRightMargin(0.0);
 594   myPadLogo->Draw();
 595   myPadLogo->cd();
 596   TASImage *myAliceLogo = new TASImage("/home/cbianchi/macros/alice_logo.png");
 597   myAliceLogo->Draw();
 598   printf("Draw logo\n");
 599   c->cd();
 600   TPaveText* t1=new TPaveText(where[0]-0.03,where[1]-0.11,where[2]+0.01,where[1]-0.01,"NDC");
 601
 602   t1->SetFillStyle(0);
 603   t1->SetBorderSize(0);
 604   //t1->AddText(0.,0.,Form("%s%s",type.Data(),(type=="Performace") ? today.Data() : ""));
 605   t1->AddText(0.,0.,Form("%s",type.Data()));
 606   t1->SetTextColor(kRed);
 607   t1->SetTextFont(42);
 608
 609   if(type=="Preliminary") {
 610     cout<<"Preliminary"<<endl;
 611     //t1->Print();
 612     t1->Draw();
 613   }
 614   else {
 615     printf("performance\n");
 616     t1->AddText(0.,0.,today.Data());
 617     t1->SetY1NDC(where[1]-0.17);
 618     t1->SetY2NDC(where[1]-0.05);
 619     t1->Draw();
 620     printf("Drawn date\n");
 621   }
 622
 623 }