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