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