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