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