7e48cc0496f9b9f4c789611c75ff34de69b5e6f0
[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 "AliEventPlaneResolutionHandler.h"
32 #include <AliVertexingHFUtils.h>
33
34 #endif
35
36 /* $Id$ */
37
38 //methods for the analysis of AliAnalysisTaskSEHFv2 output
39 //Authors: Chiara Bianchin cbianchi@pd.infn.it
40 //         Giacomo Ortona  ortona@to.infn.it
41 //         Francesco Prino prino@to.infn.it
42
43 //global variables to be set
44 TString filename="AnalysisResults_train63.root";
45 TString suffix="v2Dplus3050Cut4upcutPIDTPC";
46 TString partname="Dplus";
47 Int_t minCent=30;
48 Int_t maxCent=50;
49 //evPlane flag from AliEventPlaneResolutionHandler: 
50 //kTPCFullEta, kTPCPosEta,kVZERO,kVZEROA,kVZEROC
51 Int_t evPlane=AliEventPlaneResolutionHandler::kTPCPosEta;
52 //resolution flag fromAliEventPlaneResolutionHandler:
53 //kTwoRandSub,kTwoChargeSub,kTwoEtaSub,kThreeSub,kThreeSubTPCGap
54 Int_t evPlaneRes=AliEventPlaneResolutionHandler::kTwoEtaSub;
55 Bool_t useNcollWeight=kFALSE;
56 // pt and phi binning
57 const Int_t nptbinsnew=4;
58 Float_t ptbinsnew[nptbinsnew+1]={3.,4.,6.,8.,12.};
59 const Int_t nphibins=4;
60 Float_t phibinslim[nphibins+1]={0,TMath::Pi()/4,TMath::Pi()/2,3*TMath::Pi()/4,TMath::Pi()};
61 // mass fit configuration
62 Int_t rebin[nptbinsnew]={4,4,4,4};
63 Int_t typeb=0;  // Background: 0=expo, 1=linear, 2=pol2
64 Double_t minMassForFit=1.6;
65 Double_t maxMassForFit=2.2;
66 Float_t massRangeForCounting=0.1; // GeV
67
68 Bool_t gnopng=kFALSE; //don't save in png format (only root and eps)
69 Int_t minPtBin[nptbinsnew]={-1,-1,-1,-1};
70 Int_t maxPtBin[nptbinsnew]={-1,-1,-1,-1};
71 Double_t effInOverEffOut=1.025;
72 Double_t massD;
73
74 //methods
75 TList *LoadMassHistos(TList *inputlist,Bool_t inoutanis);
76 TList* LoadResolutionHistos(TList *inputlist);
77 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value);
78 void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs, TGraphAsymmErrors **gSignalBC1, TGraphAsymmErrors **gSignalBC2, Bool_t inoutanis);
79 void DmesonsFlowAnalysis(Bool_t inoutanis=kTRUE);
80 void DrawEventPlane();
81 Double_t GetEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3);
82
83
84
85 //method implementation
86 //_________________________________________________________________
87 Double_t GetEventPlaneResolution(Double_t& error, TH1F* hsubev1, TH1F* hsubev2, TH1F* hsubev3){
88   Double_t resolFull=1.;
89   if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoRandSub ||
90      evPlaneRes==AliEventPlaneResolutionHandler::kTwoChargeSub){
91     resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
92     error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
93   }else if(evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
94     if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta){
95       resolFull=AliVertexingHFUtils::GetFullEvResol(hsubev1);
96       error = TMath::Abs(resolFull-AliVertexingHFUtils::GetFullEvResolLowLim(hsubev1));
97     }else if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
98       resolFull=AliVertexingHFUtils::GetSubEvResol(hsubev1);
99       error = TMath::Abs(resolFull-AliVertexingHFUtils::GetSubEvResolLowLim(hsubev1));      
100     }
101   }else{
102     Double_t resolSub[3];
103     Double_t errors[3];
104     TH1F* hevplresos[3];
105     hevplresos[0]=hsubev1;
106     hevplresos[1]=hsubev2;
107     hevplresos[2]=hsubev3;
108     for(Int_t ires=0;ires<3;ires++){
109       resolSub[ires]=hevplresos[ires]->GetMean();
110       errors[ires]=hevplresos[ires]->GetMeanError();
111     }
112     Double_t lowlim[3];for(Int_t ie=0;ie<3;ie++)lowlim[ie]=TMath::Abs(resolSub[ie]-errors[ie]);
113     if(evPlane==AliEventPlaneResolutionHandler::kVZEROC ||
114        evPlane==AliEventPlaneResolutionHandler::kVZERO){
115       resolFull=TMath::Sqrt(resolSub[1]*resolSub[2]/resolSub[0]);
116       error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[1]/lowlim[0]);
117     }
118     else if(evPlane==AliEventPlaneResolutionHandler::kVZEROA){
119       resolFull=TMath::Sqrt(resolSub[0]*resolSub[2]/resolSub[1]);
120       error=resolFull-TMath::Sqrt(lowlim[2]*lowlim[0]/lowlim[1]);
121     }
122     else if(evPlane==AliEventPlaneResolutionHandler::kTPCFullEta ||
123             evPlane==AliEventPlaneResolutionHandler::kTPCPosEta){
124       resolFull=TMath::Sqrt(resolSub[0]*resolSub[1]/resolSub[2]);
125       error=resolFull-TMath::Sqrt(lowlim[0]*lowlim[1]/lowlim[2]);
126     }
127   }
128   return resolFull;
129 }
130
131 //____________________________________________________________________
132 TList* LoadResolutionHistos(TList *inputlist){
133
134   TList *outlist = new TList();
135   outlist->SetName("eventplanehistlist");
136
137   const Int_t nBins=20;
138   Double_t ncoll[nBins]={1790.77,1578.44,1394.82,1236.17,1095.08,
139                          969.86,859.571,759.959,669.648,589.588,
140                          516.039,451.409,392.853,340.493,294.426,
141                          252.385,215.484,183.284,155.101,130.963};
142
143   Int_t minCentTimesTen=minCent*10;
144   Int_t maxCentTimesTen=maxCent*10;
145   TGraphErrors* gResolVsCent=new TGraphErrors(0);
146   Int_t iPt=0;
147   Int_t nSubRes=1;
148   if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
149      evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
150   TString namereso[3]={"Reso","Reso2","Reso3"};
151   TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
152   TH2F* hevpls=(TH2F*)inputlist->FindObject(Form("hEvPlanecentr%d_%d",minCentTimesTen,minCentTimesTen+25));
153   hevpls->SetName(Form("hEvPlane%s",suffixcentr.Data()));
154   hevpls->SetTitle(Form("Event Plane angle %s",suffixcentr.Data()));
155     //    TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));QUI!!!
156   TH1F* hevplresos[3];
157   Int_t ncBin=minCentTimesTen/25;
158   
159   for(Int_t ires=0;ires<nSubRes;ires++){
160     hevplresos[ires]=(TH1F*)inputlist->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),minCentTimesTen,minCentTimesTen+25));
161     if(hevplresos[ires]){
162       hevplresos[ires]->SetName(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
163       hevplresos[ires]->SetTitle(Form("Event Plane Resolution %s%s",namereso[ires].Data(),suffixcentr.Data()));
164       if(useNcollWeight){
165         printf("Centr %d Bin %d  Ncoll %f\n",minCentTimesTen,ncBin,ncoll[ncBin]);
166         hevplresos[ires]->Scale(ncoll[ncBin]);
167       }
168     }
169   }
170   Double_t error;
171   Double_t lowestRes=1;
172   Double_t highestRes=0;
173   Double_t resolBin=GetEventPlaneResolution(error,hevplresos[0],hevplresos[1],hevplresos[2]);
174   if(resolBin<lowestRes) lowestRes=resolBin;
175   if(resolBin>highestRes) highestRes=resolBin;
176
177   Double_t binHalfWid=25./20.;
178   Double_t binCentr=(Double_t)minCentTimesTen/10.+binHalfWid;
179   gResolVsCent->SetPoint(iPt,binCentr,resolBin);
180   gResolVsCent->SetPointError(iPt,binHalfWid,error);
181   ++iPt;
182   
183   for(Int_t icentr=minCentTimesTen+25;icentr<maxCentTimesTen;icentr=icentr+25){
184     TH2F* h=(TH2F*)inputlist->FindObject(Form("hEvPlanecentr%d_%d",icentr,icentr+25));
185     if(h)hevpls->Add(h);
186     else cout<<"skipping ev plane "<<icentr<<"_"<<icentr+5<<endl;
187     TH1F* htmpresos[3];
188     for(Int_t ires=0;ires<nSubRes;ires++){
189       htmpresos[ires]=(TH1F*)inputlist->FindObject(Form("hEvPlane%scentr%d_%d",namereso[ires].Data(),icentr,icentr+25));
190       if(!htmpresos[ires])cout<<"skipping ev pl reso "<<icentr<<"_"<<icentr+25<<endl;
191     }
192     resolBin=GetEventPlaneResolution(error,htmpresos[0],htmpresos[1],htmpresos[2]);
193     if(resolBin<lowestRes) lowestRes=resolBin;
194     if(resolBin>highestRes) highestRes=resolBin;
195     binCentr=(Double_t)icentr/10.+binHalfWid;
196     gResolVsCent->SetPoint(iPt,binCentr,resolBin);
197     gResolVsCent->SetPointError(iPt,binHalfWid,error);
198     ++iPt;
199     ncBin=icentr/25;
200     for(Int_t ires=0;ires<nSubRes;ires++){
201       if(htmpresos[ires]){
202         if(useNcollWeight){
203           printf("Centr %d Bin %d  Ncoll %f\n",icentr,ncBin,ncoll[ncBin]);
204           htmpresos[ires]->Scale(ncoll[ncBin]);
205         }
206         hevplresos[ires]->Add(htmpresos[ires]);
207       }
208     }
209   }
210   outlist->Add(hevpls->Clone());
211   for(Int_t ires=0;ires<nSubRes;ires++){
212     if(hevplresos[ires]) outlist->Add(hevplresos[ires]->Clone());
213   }
214   gResolVsCent->SetName("gResolVsCent");
215   gResolVsCent->SetTitle("Resolution vs. Centrality");
216   outlist->Add(gResolVsCent->Clone());
217   return outlist;
218 }
219
220 //__________________________________________________________
221 TList *LoadMassHistos(TList *inputlist,Bool_t inoutanis){
222   // printf("Start load histos\n");
223   //  const Int_t nptbins=cutobj->GetNPtBins();
224   TList *outlist = new TList();
225   outlist->SetName("azimuthalhistoslist");
226   
227   Int_t nphi=nphibins;
228   if(inoutanis)nphi=2;
229   Int_t minCentTimesTen=minCent*10;
230   Int_t maxCentTimesTen=maxCent*10;
231   
232   //Create 2D histogram in final pt bins
233   for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
234     for(Int_t iphi=0;iphi<nphi;iphi++){
235       TH1F *hMass=0x0;//=new TH1F();
236       for(Int_t iPtBin=minPtBin[iFinalPtBin]; iPtBin<=maxPtBin[iFinalPtBin];iPtBin++){
237         for(Int_t iHisC=minCentTimesTen; iHisC<=maxCentTimesTen-25; iHisC+=25){    
238           TString hisname=Form("hMdeltaphi_pt%dcentr%d_%d",iPtBin,iHisC,iHisC+25);
239           TH2F* htmp=(TH2F*)inputlist->FindObject(hisname.Data());
240           printf("---> Histogram: %s\n",htmp->GetName());
241           Int_t startX=htmp->FindBin(phibinslim[iphi]);
242           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
243           TH1F *h1tmp;
244           if(inoutanis){
245             if(iphi==0){
246               Int_t firstBin=htmp->FindBin(0);
247               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
248               h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi0",iPtBin),firstBin,lastBin);
249               printf("In-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
250               firstBin=htmp->FindBin(3.*TMath::Pi()/4.);
251               lastBin=htmp->FindBin(TMath::Pi()-0.0001); // -0.0001 to be sure that the upper limit of the bin is pi
252               h1tmp->Add((TH1F*)htmp->ProjectionY(Form("hMass%d",iPtBin),firstBin,lastBin));
253               printf("In-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
254             }else{
255               Int_t firstBin=htmp->FindBin(TMath::Pi()/4.);
256               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
257               h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi1",iPtBin),firstBin,lastBin);
258               printf("Out-of-plane, Range: bins %d-%d -> phi %f - %f\n",firstBin,lastBin,htmp->GetXaxis()->GetBinLowEdge(firstBin),htmp->GetXaxis()->GetBinUpEdge(lastBin));
259             }
260           }else{
261             h1tmp=(TH1F*)htmp->ProjectionY(Form("hMass%d_phi%d",iPtBin,iphi),startX,endX);
262           }
263           if(hMass==0)hMass=(TH1F*)h1tmp->Clone();
264           else hMass->Add((TH1F*)h1tmp->Clone());
265         }
266       }
267       hMass->SetTitle(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
268       hMass->SetName(Form("hMass_pt%d_phi%d",iFinalPtBin,iphi));
269       outlist->Add(hMass->Clone());
270       //      hMass->DrawClone();
271       delete hMass;
272       hMass=0x0;
273     }
274   }
275   return outlist;
276 }
277 //______________________________________________________________
278 Bool_t DefinePtBins(AliRDHFCuts *cutobj){
279   Int_t nPtBinsCuts=cutobj->GetNPtBins();
280   Float_t *ptlimsCuts=cutobj->GetPtBinLimits();
281   //  for(Int_t iPt=0; iPt<nPtBinsCuts; iPt++) printf(" %d %f-%f\n",iPt,ptlimsCuts[iPt],ptlimsCuts[iPt+1]);
282   for(Int_t iPtCuts=0; iPtCuts<nPtBinsCuts; iPtCuts++){
283     for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){  
284       if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[iFinalPtBin])<0.0001){ 
285         minPtBin[iFinalPtBin]=iPtCuts;
286         if(iFinalPtBin>0) maxPtBin[iFinalPtBin-1]=iPtCuts-1;
287       }
288     }
289     if(TMath::Abs(ptlimsCuts[iPtCuts]-ptbinsnew[nptbinsnew])<0.0001) maxPtBin[nptbinsnew-1]=iPtCuts-1;
290   }
291   if(TMath::Abs(ptbinsnew[nptbinsnew]-ptlimsCuts[nPtBinsCuts])<0.0001) maxPtBin[nptbinsnew-1]=nPtBinsCuts-1;
292   for(Int_t iFinalPtBin=0; iFinalPtBin<nptbinsnew; iFinalPtBin++){
293     printf("Pt bins to be merged: %d %d\n",minPtBin[iFinalPtBin],maxPtBin[iFinalPtBin]);
294     if(minPtBin[iFinalPtBin]<0 || maxPtBin[iFinalPtBin]<0) return kFALSE;
295   }
296
297   return kTRUE;
298 }
299 //______________________________________________________________
300 Int_t GetPadNumber(Int_t ix,Int_t iy){
301   return (iy)*nptbinsnew+ix+1;
302 }
303 //______________________________________________________________
304 void FillSignalGraph(TList *histlist,TGraphAsymmErrors **gSignal,TGraphAsymmErrors **gSignalfs, TGraphAsymmErrors **gSignalBC1, TGraphAsymmErrors **gSignalBC2, Bool_t inoutanis){
305
306   Int_t nphi=nphibins;
307   if(inoutanis)nphi=2;
308   
309   //Canvases for drawing histograms
310   TCanvas *cDeltaPhi = new TCanvas("cinvmassdeltaphi","Invariant mass distributions",1200,700);
311   TCanvas *cDeltaPhifs = new TCanvas("cinvmassdeltaphifs","Invariant mass distributions - fit with fixed sigma",1200,700);
312   TCanvas *cPhiInteg = new TCanvas("cinvmass","Invariant mass distributions - #phi integrated",1200,350);
313   cDeltaPhi->Divide(nptbinsnew,nphi);
314   cDeltaPhifs->Divide(nptbinsnew,nphi);
315   cPhiInteg->Divide(nptbinsnew,1);
316   Int_t nMassBins;
317   Double_t hmin,hmax;
318   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){    
319     TH1F *histtofitfullsigma=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi0",ipt))->Clone();
320     for(Int_t iphi=0;iphi<nphi;iphi++){
321       Int_t ipad=GetPadNumber(ipt,iphi);
322       Double_t signal=0,esignal=0;
323       TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
324       if(iphi>0)histtofitfullsigma->Add((TH1F*)histtofit->Clone());
325       if(!histtofit){
326         gSignal[ipt]->SetPoint(iphi,iphi,signal);
327         gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
328         return;
329       }
330       histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
331       nMassBins=histtofit->GetNbinsX();
332       hmin=TMath::Max(minMassForFit,histtofit->GetBinLowEdge(2));
333       hmax=TMath::Min(maxMassForFit,histtofit->GetBinLowEdge(nMassBins-2));
334       histtofit->Rebin(rebin[ipt]);
335       AliHFMassFitter fitter(histtofit,hmin,hmax,1,typeb);
336       fitter.SetInitialGaussianMean(massD);
337       //      fitter.SetInitialGaussianSigma(0.012);
338       Bool_t ok=fitter.MassFitter(kFALSE);
339       if(ok){
340         fitter.DrawHere(cDeltaPhi->cd(ipad),3,1);
341         fitter.Signal(3,signal,esignal);
342       }
343       gSignal[ipt]->SetPoint(iphi,iphi,signal);
344       gSignal[ipt]->SetPointError(iphi,0,0,esignal,esignal);
345       TF1* fB1=fitter.GetBackgroundFullRangeFunc();
346       TF1* fB2=fitter.GetBackgroundRecalcFunc();
347       Float_t minBinSum=histtofit->FindBin(massD-massRangeForCounting);
348       Float_t maxBinSum=histtofit->FindBin(massD+massRangeForCounting);
349       Float_t cntSig1=0.;
350       Float_t cntSig2=0.;
351       Float_t cntErr=0.;
352       for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
353         Float_t bkg1=fB1 ? fB1->Eval(histtofit->GetBinCenter(iMB)) : 0;
354         Float_t bkg2=fB2 ? fB2->Eval(histtofit->GetBinCenter(iMB)) : 0;
355         //      printf("PtBin %d MassBin%d BC %f %f %f\n",ipt,iMB,histtofit->GetBinContent(iMB),bkg1,bkg2);
356         cntSig1+=(histtofit->GetBinContent(iMB)-bkg1);
357         cntSig2+=(histtofit->GetBinContent(iMB)-bkg2);
358         cntErr+=(histtofit->GetBinContent(iMB));
359       }
360       cntErr=TMath::Sqrt(cntErr);
361       gSignalBC1[ipt]->SetPoint(iphi,iphi,cntSig1);
362       gSignalBC1[ipt]->SetPointError(iphi,0,0,cntErr,cntErr);
363       gSignalBC2[ipt]->SetPoint(iphi,iphi,cntSig2);
364       gSignalBC2[ipt]->SetPointError(iphi,0,0,cntErr,cntErr);
365     }
366     //fit for fixed sigma
367     histtofitfullsigma->SetTitle(Form("%.1f<p_{t}<%.1f",ptbinsnew[ipt],ptbinsnew[ipt+1]));
368     nMassBins=histtofitfullsigma->GetNbinsX();
369     hmin=TMath::Max(minMassForFit,histtofitfullsigma->GetBinLowEdge(2));
370     hmax=TMath::Min(maxMassForFit,histtofitfullsigma->GetBinLowEdge(nMassBins-2));
371     histtofitfullsigma->Rebin(rebin[ipt]);
372     AliHFMassFitter fitter(histtofitfullsigma,hmin,hmax,1,typeb);
373     fitter.SetInitialGaussianMean(massD);
374     //      fitter.SetInitialGaussianSigma(0.012);
375     //fitter.SetFixGaussianSigma(0.012);
376     Bool_t ok=fitter.MassFitter(kFALSE);
377     if(ok){
378       fitter.DrawHere(cPhiInteg->cd(ipt+1),3,1);
379     }
380     Double_t sigma=fitter.GetSigma();
381     for(Int_t iphi=0;iphi<nphi;iphi++){
382       Int_t ipad=GetPadNumber(ipt,iphi);
383       TH1F *histtofit=(TH1F*)histlist->FindObject(Form("hMass_pt%d_phi%d",ipt,iphi))->Clone();
384       histtofit->SetTitle(Form("%.1f<p_{t}<%.1f, #phi%d",ptbinsnew[ipt],ptbinsnew[ipt+1],iphi));
385       nMassBins=histtofit->GetNbinsX();
386       hmin=TMath::Max(minMassForFit,histtofit->GetBinLowEdge(2));
387       hmax=TMath::Min(maxMassForFit,histtofit->GetBinLowEdge(nMassBins-2));
388       histtofit->Rebin(rebin[ipt]);
389       AliHFMassFitter fitter2(histtofit,hmin,hmax,1,typeb);
390       fitter2.SetInitialGaussianMean(massD);
391       fitter2.SetFixGaussianSigma(sigma);
392       Bool_t ok2=fitter2.MassFitter(kFALSE);
393       Double_t signal=0,esignal=0;
394       if(ok2){
395         fitter2.DrawHere(cDeltaPhifs->cd(ipad),3,1);
396         fitter2.Signal(3,signal,esignal);
397       }
398       gSignalfs[ipt]->SetPoint(iphi,iphi,signal);
399       gSignalfs[ipt]->SetPointError(iphi,0,0,esignal,esignal);
400     }
401   }//end loop on pt bin
402
403
404   cDeltaPhi->SaveAs(Form("InvMassDeltaPhi_%s.eps",suffix.Data()));
405   cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.eps",suffix.Data()));
406   cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.root",suffix.Data()));
407   cPhiInteg->SaveAs(Form("InvMassfullphi_%s.eps",suffix.Data()));
408
409   if(!gnopng){
410     cDeltaPhi->SaveAs(Form("InvMassDeltaPhi_%s.png",suffix.Data()));
411     cDeltaPhifs->SaveAs(Form("InvMassDeltaPhi_fs_%s.png",suffix.Data()));
412     cPhiInteg->SaveAs(Form("InvMassfullphi_%s.png",suffix.Data()));
413   }
414   //cDeltaPhifs->DrawClone();
415   //cDeltaPhifs->Close();
416  
417 }
418 //______________________________________________________________
419 void DmesonsFlowAnalysis(Bool_t inoutanis){
420   TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
421   TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
422
423   AliRDHFCuts *cutsobj=0x0;
424   //Load input data from AliAnalysisTaskSEHFv2
425   TFile *f = TFile::Open(filename.Data());
426   if(!f){
427     printf("file %s not found, please check file name\n",filename.Data());return;
428   }
429   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
430   if(!dir){
431     printf("Directory %s not found, please check dir name\n",dirname.Data());return;
432   }
433   if(partname.Contains("Dzero")) {
434     cutsobj=((AliRDHFCutsD0toKpi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
435     massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
436   }
437   if(partname.Contains("Dplus")){
438     cutsobj=((AliRDHFCutsDplustoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
439     massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
440   }
441   if(partname.Contains("Dstar")) {
442     cutsobj=((AliRDHFCutsDStartoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
443     massD=(TDatabasePDG::Instance()->GetParticle(413)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass()); 
444   }
445
446   TList *list =(TList*)dir->Get(listname.Data());
447   if(!list){
448     printf("list %s not found in file, please check list name\n",listname.Data());return;
449   }
450   if(!cutsobj){
451     printf("cut object not found in file, please check keylist number\n");return;
452   }
453   //Define new pt bins
454   if(!DefinePtBins(cutsobj)){
455     printf("cut not define pt bins\n");return;
456   }
457   
458   //Load mass histograms corresponding to the required centrality, pt range and phi binning
459   printf("Load mass histos \n");
460   TList *histlist=LoadMassHistos(list,inoutanis);
461   TString aniss="";
462   if(inoutanis)aniss+="anis";
463   histlist->SaveAs(Form("v2Histograms_%d_%d_%s_%s.root",minCent,maxCent,aniss.Data(),suffix.Data()),"RECREATE");
464
465   Int_t nphi=nphibins;
466   if(inoutanis)nphi=2;
467
468
469   printf("average pt for pt bin \n");
470   //average pt for pt bin
471   AliVertexingHFUtils *utils=new AliVertexingHFUtils();
472   Int_t minCentTimesTen=minCent*10;
473   Int_t maxCentTimesTen=maxCent*10;
474   TH2F* hmasspt=(TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",minCentTimesTen,minCentTimesTen+25));
475   for(Int_t icent=minCentTimesTen+25;icent<maxCentTimesTen;icent=icent+25)hmasspt->Add((TH2F*)list->FindObject(Form("hMPtCandcentr%d_%d",icent,icent+25)));
476   Float_t averagePt[nptbinsnew];
477   Float_t errorPt[nptbinsnew];
478   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
479     Int_t binMin=hmasspt->FindBin(ptbinsnew[ipt]);
480     Int_t binMax=hmasspt->FindBin(ptbinsnew[ipt+1]-0.001);
481     if(TMath::Abs(hmasspt->GetXaxis()->GetBinLowEdge(binMin)-ptbinsnew[ipt])>0.001 || 
482        TMath::Abs(hmasspt->GetXaxis()->GetBinUpEdge(binMax)-ptbinsnew[ipt+1])>0.001){
483       printf("Error in pt bin limits for projection!\n");
484       return;
485     }
486     TH1F *histtofit = (TH1F*)hmasspt->ProjectionY("_py",binMin,binMax);
487     Int_t nMassBins=histtofit->GetNbinsX();
488     Double_t hmin=histtofit->GetBinLowEdge(2); // need wide range for <pt>
489     Double_t hmax=histtofit->GetBinLowEdge(nMassBins-2); // need wide range for <pt>
490     AliHFMassFitter fitter(histtofit,hmin,hmax,1);
491     fitter.MassFitter(kFALSE);
492     Float_t massFromFit=fitter.GetMean();
493     Float_t sigmaFromFit=fitter.GetSigma();
494     TF1* funcB2=fitter.GetBackgroundRecalcFunc();
495     utils->AveragePt(averagePt[ipt],errorPt[ipt],ptbinsnew[ipt],ptbinsnew[ipt+1],hmasspt,massFromFit,sigmaFromFit,funcB2,2.5,4.5,1);
496   }
497   printf("Average pt\n");
498   for(Int_t ipt=0;ipt<nptbinsnew;ipt++) printf("%f +- %f\n",averagePt[ipt],errorPt[ipt]);
499
500   printf("Fill TGraphs for signal \n");
501   //Fill TGraphs for signal
502   TGraphAsymmErrors *gSignal[nptbinsnew];
503   TGraphAsymmErrors *gSignalfs[nptbinsnew];
504   TGraphAsymmErrors *gSignalBC1[nptbinsnew];
505   TGraphAsymmErrors *gSignalBC2[nptbinsnew];
506   for(Int_t i=0;i<nptbinsnew;i++){
507     gSignal[i]=new TGraphAsymmErrors(nphi);
508     gSignal[i]->SetName(Form("gasigpt%d",i));
509     gSignal[i]->SetTitle(Form("Signal %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
510     gSignal[i]->SetMarkerStyle(25);
511     gSignalfs[i]=new TGraphAsymmErrors(nphi);
512     gSignalfs[i]->SetName(Form("gasigfspt%d",i));
513     gSignalfs[i]->SetTitle(Form("Signal (fixed sigma) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
514     gSignalfs[i]->SetMarkerStyle(21);
515     gSignalBC1[i]=new TGraphAsymmErrors(nphi);
516     gSignalBC1[i]->SetName(Form("gasigBC1pt%d",i));
517     gSignalBC1[i]->SetTitle(Form("Signal (BC1) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
518     gSignalBC2[i]=new TGraphAsymmErrors(nphi);
519     gSignalBC2[i]->SetName(Form("gasigBC2pt%d",i));
520     gSignalBC2[i]->SetTitle(Form("Signal (BC2) %.1f<p_{t}<%.1f GeV/c;#Delta#phi bin;Counts",ptbinsnew[i],ptbinsnew[i+1]));
521   }
522   FillSignalGraph(histlist,gSignal,gSignalfs,gSignalBC1,gSignalBC2,inoutanis);
523
524   //EP resolution
525   // TList* resolhist=LoadResolutionHistos(list);
526   // TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
527   // TH1F* hevplresos[3];
528   // TString namereso[3]={"Reso","Reso2","Reso3"};
529   // Int_t nSubRes=1;
530   // if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
531   //    evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
532   // for(Int_t ires=0;ires<nSubRes;ires++){
533   //   hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
534   // }
535   // Double_t errorres;
536   // Double_t resol=GetEventPlaneResolution(errorres,hevplresos[0],hevplresos[1],hevplresos[2]);
537   AliEventPlaneResolutionHandler* epres=new AliEventPlaneResolutionHandler();
538   epres->SetEventPlane(evPlane);
539   epres->SetResolutionOption(evPlaneRes);
540   epres->SetUseNcollWeights(useNcollWeight);
541   Double_t resol=epres->GetEventPlaneResolution(minCent,maxCent);
542   delete epres;
543   printf("Event plane resolution %f\n",resol);
544
545   printf("Compute v2 \n");
546   //compute v2
547   TCanvas *cv2fs =new TCanvas("v2_fs","v2 Fit with fixed sigma");
548   TCanvas *cv2 =new TCanvas("v2","v2 - systematic on yield extraction");
549   TGraphAsymmErrors *gv2=new TGraphAsymmErrors(nptbinsnew);
550   TGraphAsymmErrors *gv2fs=new TGraphAsymmErrors(nptbinsnew);
551   TGraphAsymmErrors *gv2BC1=new TGraphAsymmErrors(nptbinsnew);
552   TGraphAsymmErrors *gv2BC2=new TGraphAsymmErrors(nptbinsnew);
553   TGraphAsymmErrors *grelSystEff=new TGraphAsymmErrors(nptbinsnew);
554
555   gv2->SetName("gav2");
556   gv2->SetTitle("v_{2}(p_{t})");
557   gv2fs->SetName("gav2fs");
558   gv2fs->SetTitle("v_{2}(p_{t}) (fit with fixed sigma)");
559   gv2BC1->SetName("gav2BC1");
560   gv2BC1->SetTitle("v_{2}(p_{t}) (bin counting 1)");
561   gv2BC2->SetName("gav2BC2");
562   gv2BC2->SetTitle("v_{2}(p_{t}) (bin counting 2)");
563   grelSystEff->SetName("grelSystEff");
564   grelSystEff->SetTitle("SystErrEff");
565
566   //Prepare output file
567   TFile *fout=new TFile(Form("v2Output_%d_%d_%s_%s.root",minCent,maxCent,aniss.Data(),suffix.Data()),"RECREATE");
568
569   for(Int_t ipt=0;ipt<nptbinsnew;ipt++){
570     fout->cd();
571     gSignal[ipt]->Write();
572     gSignalfs[ipt]->Write();
573     gSignalBC1[ipt]->Write();
574     gSignalBC2[ipt]->Write();
575
576     if(inoutanis){
577       //v2 from in-out anisotropy
578       Double_t *y,*yfs,*ybc1,*ybc2;
579       y=gSignal[ipt]->GetY();
580       yfs=gSignalfs[ipt]->GetY();
581       ybc1=gSignalBC1[ipt]->GetY();
582       ybc2=gSignalBC2[ipt]->GetY();
583
584       Double_t nIn=y[0];
585       Double_t nOut=y[1];
586       Double_t enIn=gSignal[ipt]->GetErrorY(0);
587       Double_t enOut=gSignal[ipt]->GetErrorY(1);
588       Double_t anis=(nIn-nOut)/(nIn+nOut);
589       //      Double_t eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
590       Double_t eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
591       Double_t v2=anis*TMath::Pi()/4./resol;
592       Double_t ev2=eAnis*TMath::Pi()/4./resol;
593       gv2->SetPoint(ipt,averagePt[ipt],v2);
594       gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
595
596       nIn=yfs[0];
597       nOut=yfs[1];
598       enIn=gSignalfs[ipt]->GetErrorY(0);
599       enOut=gSignalfs[ipt]->GetErrorY(1);
600       anis=(nIn-nOut)/(nIn+nOut);
601       //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
602       eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
603       v2=anis*TMath::Pi()/4./resol;
604       ev2=eAnis*TMath::Pi()/4./resol;
605       gv2fs->SetPoint(ipt,averagePt[ipt],v2);
606       gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
607       Double_t anis1=(nIn-nOut*effInOverEffOut)/(nIn+nOut*effInOverEffOut);
608       Double_t anis2=(nIn*effInOverEffOut-nOut)/(nIn*effInOverEffOut+nOut);
609       Double_t systEffUp=0.,systEffDown=0.;
610       if(anis1>anis && anis1>anis2) systEffUp=anis1/anis;
611       if(anis2>anis && anis2>anis1) systEffUp=anis2/anis;
612       if(anis1<anis && anis1<anis2) systEffDown=anis1/anis;
613       if(anis2<anis && anis2<anis1) systEffDown=anis2/anis;
614       printf(" Bin %d <pt>=%.3f  v2=%f+-%f systEff=%f %f\n",ipt,averagePt[ipt],v2,ev2,systEffUp*v2,systEffDown*v2);
615       grelSystEff->SetPoint(ipt,averagePt[ipt],v2);
616       grelSystEff->SetPointError(ipt,0.4,0.4,v2*(1-systEffDown),v2*(systEffUp-1));
617
618       nIn=ybc1[0];
619       nOut=ybc1[1];
620       enIn=gSignalBC1[ipt]->GetErrorY(0);
621       enOut=gSignalBC1[ipt]->GetErrorY(1);
622       anis=(nIn-nOut)/(nIn+nOut);
623       //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
624       eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
625       v2=anis*TMath::Pi()/4./resol;
626       ev2=eAnis*TMath::Pi()/4./resol;
627       gv2BC1->SetPoint(ipt,averagePt[ipt],v2);
628       gv2BC1->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
629
630       nIn=ybc2[0];
631       nOut=ybc2[1];
632       enIn=gSignalBC2[ipt]->GetErrorY(0);
633       enOut=gSignalBC2[ipt]->GetErrorY(1);
634       anis=(nIn-nOut)/(nIn+nOut);
635       //     eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
636       eAnis=2./((nIn+nOut)*(nIn+nOut))*TMath::Sqrt(nIn*nIn*enOut*enOut+nOut*nOut*enIn*enIn);
637       v2=anis*TMath::Pi()/4./resol;
638       ev2=eAnis*TMath::Pi()/4./resol;
639       gv2BC2->SetPoint(ipt,averagePt[ipt],v2);
640       gv2BC2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
641
642    }else{
643       TF1 *flowFunc = new TF1("flow","[0]*(1.+2.*[1]*TMath::Cos(2.*x))");
644       //v2 from fit to Deltaphi distribution
645       gSignal[ipt]->Fit(flowFunc);
646       Double_t v2 = flowFunc->GetParameter(1)/resol;
647       Double_t ev2=flowFunc->GetParError(1)/resol;
648       gv2->SetPoint(ipt,averagePt[ipt],v2);
649       gv2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
650
651       gSignalfs[ipt]->Fit(flowFunc);
652       v2 = flowFunc->GetParameter(1)/resol;
653       ev2=flowFunc->GetParError(1)/resol;
654       gv2fs->SetPoint(ipt,averagePt[ipt],v2);
655       gv2fs->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
656
657       gSignalBC1[ipt]->Fit(flowFunc);
658       v2 = flowFunc->GetParameter(1)/resol;
659       ev2=flowFunc->GetParError(1)/resol;
660       gv2BC1->SetPoint(ipt,averagePt[ipt],v2);
661       gv2BC1->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
662
663       gSignalBC2[ipt]->Fit(flowFunc);
664       v2 = flowFunc->GetParameter(1)/resol;
665       ev2=flowFunc->GetParError(1)/resol;
666       gv2BC2->SetPoint(ipt,averagePt[ipt],v2);
667       gv2BC2->SetPointError(ipt,averagePt[ipt]-ptbinsnew[ipt],ptbinsnew[ipt+1]-averagePt[ipt],ev2,ev2);
668     }
669   }
670
671   gv2->SetMarkerStyle(22);
672   gv2->SetMarkerColor(2);
673   gv2->SetLineColor(2);
674   gv2fs->SetMarkerStyle(20);
675   gv2fs->SetMarkerColor(1);
676   gv2fs->SetLineColor(1);
677   gv2BC1->SetMarkerStyle(23);
678   gv2BC1->SetMarkerColor(kGreen+1);
679   gv2BC1->SetLineColor(kGreen+1);
680   gv2BC2->SetMarkerStyle(26);
681   gv2BC2->SetMarkerColor(kGreen+1);
682   gv2BC2->SetLineColor(kGreen+1);
683   
684   cv2->cd();
685   gv2fs->Draw("AP");
686   gv2fs->SetMinimum(-0.2);
687   gv2fs->SetMaximum(0.6);
688   gv2fs->GetYaxis()->SetTitle("v_{2}");
689   gv2fs->GetXaxis()->SetTitle("p_{t} [GeV/c]");
690   gv2->Draw("PSAME");
691   gv2BC1->Draw("PSAME");
692   gv2BC2->Draw("PSAME");
693   TLegend* leg=new TLegend(0.65,0.7,0.89,0.89);
694   leg->SetFillStyle(0);
695   leg->SetBorderSize(0);
696   leg->SetTextFont(42);
697   leg->AddEntry(gv2fs,"Fixed sigma","P")->SetTextColor(gv2fs->GetMarkerColor());
698   leg->AddEntry(gv2,"Free sigma","P")->SetTextColor(gv2->GetMarkerColor());
699   leg->AddEntry(gv2BC1,"Bin counting 1","P")->SetTextColor(gv2BC1->GetMarkerColor());
700   leg->AddEntry(gv2BC2,"Bin counting 2","P")->SetTextColor(gv2BC2->GetMarkerColor());
701   leg->Draw();
702   cv2fs->cd();
703   gv2fs->Draw("AP");
704   gv2fs->GetYaxis()->SetTitle("v_{2}");
705   gv2fs->GetXaxis()->SetTitle("p_{t} [GeV/c]");
706
707   fout->cd();
708   gv2->Write();
709   gv2fs->Write();
710   gv2BC1->Write();
711   gv2BC2->Write();
712   grelSystEff->Write();
713 }
714 //___________________________________________________________
715 Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value){
716   for (Int_t i=0;i<nbins;i++){
717     if(value>=array[i] && value<array[i+1]){
718       return i;
719     }
720   }
721   cout<<value<< " out of range "<<array[0]<<", "<<array[nbins]<<endl;
722   return -1;
723 }
724
725 //__________________________________________________________
726 void DrawEventPlane(){
727   TString dirname=Form("PWGHF_D2H_HFv2_%s%s",partname.Data(),suffix.Data());
728   TString listname=Form("coutputv2%s%s",partname.Data(),suffix.Data());
729   TFile *f = TFile::Open(filename.Data());
730   if(!f){
731     printf("file %s not found, please check file name\n",filename.Data());return;
732   }
733   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname.Data());
734   if(!dir){
735     printf("Directory %s not found, please check dir name\n",dirname.Data());return;
736   }
737   TList *list =(TList*)dir->Get(listname.Data());
738   if(!list){
739     printf("list %s not found in file, please check list name\n",listname.Data());return;
740   } 
741   if(evPlane==AliEventPlaneResolutionHandler::kTPCPosEta &&
742      evPlaneRes==AliEventPlaneResolutionHandler::kTwoEtaSub){
743     printf("Wrong setting of event plane resolution forced it to kTwoSub\n");
744     evPlaneRes=AliEventPlaneResolutionHandler::kTwoRandSub;
745   }
746   TList* resolhist=LoadResolutionHistos(list);
747   TGraphErrors* gResolVsCent=(TGraphErrors*)resolhist->FindObject("gResolVsCent");
748   Double_t lowestRes=1;
749   Double_t highestRes=0;
750   for(Int_t ipt=0; ipt<gResolVsCent->GetN(); ipt++){    
751     Double_t x,resolBin;
752     gResolVsCent->GetPoint(ipt,x,resolBin);
753     if(resolBin<lowestRes) lowestRes=resolBin;
754     if(resolBin>highestRes) highestRes=resolBin;
755   }
756   TString suffixcentr=Form("centr%d_%d",minCent,maxCent);
757   TH2F* hevpls=(TH2F*)resolhist->FindObject(Form("hEvPlane%s",suffixcentr.Data()));
758   TH1F* hevplresos[3];
759   TString namereso[3]={"Reso","Reso2","Reso3"};
760   Int_t nSubRes=1;
761   if(evPlaneRes==AliEventPlaneResolutionHandler::kThreeSub||
762      evPlaneRes==AliEventPlaneResolutionHandler::kThreeSubTPCGap)nSubRes=3;
763   for(Int_t ires=0;ires<nSubRes;ires++){
764     hevplresos[ires]=(TH1F*)resolhist->FindObject(Form("hEvPlane%s%s",namereso[ires].Data(),suffixcentr.Data()));
765   }
766
767   gStyle->SetCanvasColor(0);
768   gStyle->SetTitleFillColor(0);
769   gStyle->SetStatColor(0);
770   gStyle->SetOptStat(0);
771   //gStyle->SetPalette(1);
772   gStyle->SetOptFit(1);
773   gStyle->SetLabelFont(42,"XYZ");
774   gStyle->SetTextFont(42);
775   gStyle->SetStatFont(42);
776   //  gStyle->SetTitleFont(42);
777   gStyle->SetStatY(0.89);
778   gStyle->SetStatX(0.89);
779   //  gStyle->SetTitleAlign(22);
780   gStyle->SetTitleFont(42,"xyzg");
781
782   TH1F* htpc = (TH1F*)hevpls->ProjectionX();
783   TH1F* hv0 = (TH1F*)hevpls->ProjectionY();
784   TH1F* hUsedPlane;
785   if(evPlane==AliEventPlaneResolutionHandler::kVZERO ||
786      evPlane==AliEventPlaneResolutionHandler::kVZEROA ||
787      evPlane==AliEventPlaneResolutionHandler::kVZEROC) hUsedPlane=(TH1F*)hv0->Clone("hUsedPlane");
788   else hUsedPlane=(TH1F*)htpc->Clone("hUsedPlane");
789
790   TCanvas* cvtotevpl=new TCanvas(Form("cvtotevpl%s",suffixcentr.Data()),Form("Ev plane %s",suffixcentr.Data()),1200,400);
791   cvtotevpl->Divide(3,1);
792   cvtotevpl->cd(1);
793   hevpls->SetTitleFont(42);
794   hevpls->Draw("COLZ");
795   cvtotevpl->cd(2);
796   htpc->Draw();
797   htpc->Fit("pol0");
798   cvtotevpl->cd(3);
799   hv0->Draw();
800   hv0->Fit("pol0");
801   
802   TCanvas* cvfitevpl=new TCanvas(Form("cvditevpl%s",suffixcentr.Data()),Form("Fit to Ev plane %s",suffixcentr.Data()));
803   cvfitevpl->cd();
804   hUsedPlane->Draw();
805   TF1* four2=new TF1("four2","[0]*(1+2*[1]*cos(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
806   TF1* four2s=new TF1("four2s","[0]*(1+2*[1]*sin(2*x))",hUsedPlane->GetXaxis()->GetXmin(),hUsedPlane->GetXaxis()->GetXmax());
807   hUsedPlane->Fit(four2s);
808   hUsedPlane->Fit(four2);
809   four2->SetLineColor(1);
810   four2s->SetLineColor(4);
811   four2->Draw("same");
812   four2s->Draw("same");
813   hUsedPlane->SetMaximum(hUsedPlane->GetMaximum()*1.07);
814   TLatex* tsin=new TLatex(0.15,0.84,Form("1+2*(%.4f)*sin(2*#Phi_{EP})",four2s->GetParameter(1)));
815   tsin->SetNDC();
816   tsin->SetTextColor(4);
817   tsin->Draw();
818   TLatex* tcos=new TLatex(0.15,0.77,Form("1+2*(%.4f)*cos(2*#Phi_{EP})",four2->GetParameter(1)));
819   tcos->SetNDC();
820   tcos->SetTextColor(1);
821   tcos->Draw();
822   
823   // Compute the second Fourier component for since and cosine
824   Double_t aveCos2Phi=0.;
825   Double_t aveSin2Phi=0.;
826   Double_t counts=0.;
827   for(Int_t i=1; i<=hUsedPlane->GetNbinsX(); i++){
828     Double_t co=hUsedPlane->GetBinContent(i);
829     Double_t phi=hUsedPlane->GetBinCenter(i);
830     aveCos2Phi+=co*TMath::Cos(2*phi);
831     aveSin2Phi+=co*TMath::Sin(2*phi);
832     counts+=co;
833   }
834   if(counts>0){ 
835     aveCos2Phi/=counts;
836     aveSin2Phi/=counts;
837   }
838   printf("\n------ Fourier 2nd components of EP distribution ------\n");
839   printf("<cos(2*Psi_EP)>=%f\n",aveCos2Phi);
840   printf("<sin(2*Psi_EP)>=%f\n",aveSin2Phi);
841   printf("\n");
842
843
844   TCanvas* cvtotevplreso=new TCanvas(Form("cvtotevplreso%s",suffixcentr.Data()),Form("Ev plane Resolution %s",suffixcentr.Data()));
845   cvtotevplreso->cd();
846   hevplresos[0]->SetTitle("TPC cos2(#Psi_{A}-#Psi_{B})");
847   hevplresos[0]->Draw();
848   gStyle->SetOptStat(0);
849   if(nSubRes>1){
850     hevplresos[1]->SetLineColor(2);
851     hevplresos[2]->SetLineColor(3);
852     hevplresos[0]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0A})");
853     hevplresos[1]->SetTitle("cos2(#Psi_{TPC}-#Psi_{V0C})");
854     hevplresos[2]->SetTitle("cos2(#Psi_{V0A}-#Psi_{V0C})");
855     hevplresos[1]->Draw("SAME");
856     hevplresos[2]->Draw("SAME");
857   }
858   TLegend *leg = cvtotevplreso->BuildLegend();
859   leg->SetLineColor(0);
860   leg->SetFillColor(0);
861
862   Double_t error;
863   Double_t resolFull=GetEventPlaneResolution(error,hevplresos[0],hevplresos[1],hevplresos[2]);
864   
865   TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
866   pvreso->SetBorderSize(0);
867   pvreso->SetFillStyle(0);
868   pvreso->AddText(Form("Resolution on full event = %.4f#pm%.4f\n",resolFull,error));
869   pvreso->Draw();
870
871   gResolVsCent->SetMarkerStyle(20);
872   gResolVsCent->GetXaxis()->SetTitle("Centrality");
873   gResolVsCent->GetYaxis()->SetTitle("EP Resolution");
874  
875   Int_t colorave=kBlue-7;
876   if(useNcollWeight) colorave=kOrange+1;
877
878   TCanvas* cresvscent=new TCanvas("cresvscent","ResolVsCent");
879   cresvscent->cd();
880   gResolVsCent->Draw("AP");
881   TLine* lin=new TLine(gResolVsCent->GetXaxis()->GetXmin(),resolFull,gResolVsCent->GetXaxis()->GetXmax(),resolFull);
882   lin->SetLineColor(colorave);
883   lin->SetLineWidth(2);
884   lin->Draw();
885   TLatex* taveres=new TLatex(gResolVsCent->GetXaxis()->GetXmin()+0.5,resolFull*1.01,Form("<R_{2}>=%.4f",resolFull));
886   taveres->SetTextColor(colorave);
887   taveres->Draw();
888   printf("\n\n---- Syst from centrality dependence of resolution ----\n");
889   printf("MinRes=%f  MaxRes=%f AveRes=%f ---> Syst: -%f%% +%f%%\n",lowestRes,highestRes,resolFull,(resolFull-lowestRes)/resolFull,(highestRes-resolFull)/resolFull);
890
891   TFile* fout=new TFile(Form("EvPlanecentr%d-%d.root",minCent,maxCent),"recreate");
892   fout->cd();
893   hevpls->Write();
894   for(Int_t ires=0;ires<nSubRes;ires++)hevplresos[ires]->Write();
895   gResolVsCent->Write();
896
897
898 }
899