From Francesco: Added AliAnalysisTaskVnV0 and macros
[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / extractFlowVZERO.C
1 extractFlowVZERO(Int_t itech=0,Int_t ibin = 3,Int_t step = 0,Int_t arm=2,Float_t pTh = 0.9,Int_t charge=0,Int_t addbins=0){
2   TCanvas *can1 = new TCanvas("cV0Acheck","cV0Acheck");
3
4 /*
5 read outVZEROv2.root or outVZEROv3.root files
6
7 itech=0 TPC stand alone PID and selection
8 itech=1 TPC stand alone PID but TOF track selection
9 itech=2 TPC&TOF PID (TOF strictly required)
10 itech=3 TPC|TOF PID (the same of before but using TPC stand alone when TOF not available)
11
12 ibin = centrality bin
13
14 step = species (0=all charges, 1=pi, 2=K, 3=p, 4=e, 5=d, 6=t, 7=3He
15
16 arm = armonic 2 (v2) or 3 (v3)
17
18 pTh = probability threshold = 0.6 , 0.8 , 0.9 , 0.95
19
20 charge = 1(pos), -1(neg), 0(both)
21
22 addbins = to merge more centrality bin (i.e. ibin = 2, addbins = 3 merge 2,3,4,5 centrality bins) - preliminary version!!!!!!
23
24 */
25
26   // load libraries
27   gSystem->Load("libANALYSIS");
28   gSystem->Load("libANALYSISalice");
29   gSystem->Load("libCORRFW.so");
30
31   // get objects
32   char filein[100];
33   sprintf(filein,"outVZEROv%i.root",arm);
34   TFile *f=new TFile(filein);
35   sprintf(filein,"contVZEROv%i",arm);
36   TList *l = f->Get(filein);
37   l->ls();
38   AliCFContainer *c = l->At(0);
39   AliCFContainer *c3 = l->At(1);
40   TProfile *p1 = l->At(2);
41   TProfile *p2 = l->At(3);
42   TProfile *p3 = l->At(4);
43   Int_t iVar[2] = {1,3};
44
45   Int_t iVarPsi[1] = {5};
46   /* 0=centrality , 1= cos(arm*deltaphi) , 2=charge , 3=pt , 4=probability , 5=EP , 6=PIDmask*/
47   Double_t iMin[7] = {ibin,-1,-1.5,0,pTh+0.00001,-3.14,0};
48   Double_t iMax[7] = {ibin+addbins,1,1.5,20,1.1,3.14,1};
49
50   char *iTechName[4] = {"tpc","tpcInTof","tof","tpcORtof"};
51
52   Bool_t kOnlyTPC = kFALSE;
53
54   if(itech==0 && (step!=0)){
55       kOnlyTPC = kTRUE;
56   }
57   else if(itech==1 && (step!=0)){
58       iMin[6] = 2;
59       iMax[6] = 2;
60   }
61   else if(itech==2 || itech==1){
62       iMin[6] = 1;
63       iMax[6] = 1;
64   }
65   else{
66       iMin[6] = 0;
67       iMax[6] = 1;
68   }
69
70   if(charge==-1) iMax[2]=0;
71   if(charge==1) iMin[2]=0;
72
73   // EP distribution needed for NUA corrections
74   TH2F *hPhiA = l->At(5);
75   TH2F *hPhiC = l->At(6);
76
77   TH1D *hPhiPA = hPhiA->ProjectionY("_py",iMin[0]+1,iMax[0]+1);
78   TH1D *hPhiPC = hPhiC->ProjectionY("_py",iMin[0]+1,iMax[0]+1);
79
80   // EP resolution variables
81   Float_t res1=0,res2=0,res3=0; 
82   Float_t eres1=0,eres2=0,eres3=0; 
83
84   for(Int_t i=iMin[0];i<=iMax[0];i++){ // in case of more centrality bins weighting with the errors
85     if(p1->GetBinError(i+1)){
86       eres1 += 1./p1->GetBinError(i+1)/p1->GetBinError(i+1);
87       res1 += p1->GetBinContent(i+1)/p1->GetBinError(i+1)/p1->GetBinError(i+1);      
88     }
89     if(p2->GetBinError(i+1)){
90       eres2 += 1./p2->GetBinError(i+1)/p2->GetBinError(i+1);
91       res2 += p2->GetBinContent(i+1)/p2->GetBinError(i+1)/p2->GetBinError(i+1);      
92     }
93     if(p3->GetBinError(i+1)){
94       eres3 += 1./p3->GetBinError(i+1)/p3->GetBinError(i+1);
95       res3 += p3->GetBinContent(i+1)/p3->GetBinError(i+1)/p3->GetBinError(i+1);      
96     }
97
98     if(eres1) res1 /= eres1;
99     if(eres2) res2 /= eres2;
100     if(eres3) res3 /= eres3;
101
102     if(eres1) eres1 = TMath::Sqrt(1./eres1);
103     if(eres1) eres2 = TMath::Sqrt(1./eres2);
104     if(eres1) eres3 = TMath::Sqrt(1./eres3);
105   }
106
107   // NUA correction (fit to EP distribution)
108   TF1 *fpol = new TF1("fPol","pol0");
109   hPhiPA->Fit("fPol","","",-TMath::Pi()/arm,TMath::Pi()/arm);
110   Float_t scalA = fPol->GetParameter(0);
111   hPhiPC->Fit("fPol","","",-TMath::Pi()/arm,TMath::Pi()/arm);
112   Float_t scalC = fPol->GetParameter(0);
113
114   AliCFContainer *c2[20];
115   AliCFContainer *c4[20];
116   TH2D *h[20],*h2[20];
117
118   AliCFContainer *c2bis[20];
119   AliCFContainer *c4bis[20];
120
121   AliCFContainer *cPsi2[20];
122   AliCFContainer *cPsi4[20];
123   TH1D *hPsi[20],*hPsi2[20];
124   AliCFContainer *cPsi2bis[20];
125   AliCFContainer *cPsi4bis[20];
126
127   Float_t intA = 0;
128   Float_t intC = 0;
129
130   for(Int_t i=5;i<15;i++){
131     printf("%i\n",i);
132     iMin[5] = hPhiPA->GetBinCenter(i+1);
133     iMax[5] = hPhiPA->GetBinCenter(i+1);
134     if(!kOnlyTPC){
135       c2[i] = c->MakeSlice(2,iVar,iMin,iMax);
136       c4[i]= c3->MakeSlice(2,iVar,iMin,iMax);
137     }
138     else{ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF)
139       iMin[6] = 0;
140       iMax[6] = 0;
141       c2[i] = c->MakeSlice(2,iVar,iMin,iMax);
142       c4[i]= c3->MakeSlice(2,iVar,iMin,iMax);
143
144       iMin[6] = 2;
145       iMax[6] = 2;
146       c2bis[i] = c->MakeSlice(2,iVar,iMin,iMax);
147       c4bis[i]= c3->MakeSlice(2,iVar,iMin,iMax);
148    }
149
150
151     if(!kOnlyTPC){
152       cPsi2[i] = c->MakeSlice(1,iVarPsi,iMin,iMax);
153       cPsi4[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax);
154     }
155     else{ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF)
156       iMin[6] = 0;
157       iMax[6] = 0;
158       cPsi2[i] = c->MakeSlice(1,iVarPsi,iMin,iMax);
159       cPsi4[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax);
160
161       iMin[6] = 2;
162       iMax[6] = 2;
163       cPsi2bis[i] = c->MakeSlice(1,iVarPsi,iMin,iMax);
164       cPsi4bis[i]= c3->MakeSlice(1,iVarPsi,iMin,iMax);
165     }
166
167     h[i] = (TH2D *) c2[i]->Project(step,0,1);
168     h2[i] = (TH2D *) c4[i]->Project(step,0,1);
169     if(kOnlyTPC){
170       TH2D *htemp = (TH2D *) c2bis[i]->Project(step,0,1);
171       h[i]->Add(htemp);
172       htemp = (TH2D *) c4bis[i]->Project(step,0,1);
173       h2[i]->Add(htemp);
174     }
175
176     if(hPhiPA->GetBinContent(i)) h[i]->Scale(scalA/hPhiPA->GetBinContent(i+1)); // reweighting for NUA correction
177     if(hPhiPC->GetBinContent(i)) h2[i]->Scale(scalC/hPhiPC->GetBinContent(i+1)); // reweighting for NUA correction
178
179     hPsi[i] = (TH1D *) cPsi2[i]->Project(step,0);
180     hPsi2[i] = (TH1D *) cPsi4[i]->Project(step,0);
181
182     if(kOnlyTPC){ // merge to maskPID bins (TPC standalone without TOF + TPC stand alone with TOF)
183       TH1D *htemp2 = (TH1D *) cPsi2bis[i]->Project(step,0);
184       hPsi[i]->Add(htemp2);
185       htemp2 = (TH1D *) cPsi4bis[i]->Project(step,0);
186       hPsi2[i]->Add(htemp2);
187     }
188
189     if(hPhiPA->GetBinContent(i+1)) hPsi[i]->Scale(scalA/hPhiPA->GetBinContent(i+1)); // check of reweighting for NUA correction
190     if(hPhiPC->GetBinContent(i+1)) hPsi2[i]->Scale(scalC/hPhiPC->GetBinContent(i+1)); // check of reweighting for NUA correction
191
192     intA += hPsi[i]->Integral();
193     intC += hPsi2[i]->Integral();
194
195   }
196
197 //   hPhiPA->Scale(intA / (scalA)/10);
198 //   hPhiPC->Scale(intC / (scalC)/10);
199
200   // NUA correction check V0A
201   hPhiPA->Draw();
202   hPsi[5]->Draw("SAME");
203   hPsi[5]->Scale(scalA / intA * 10);
204   for(Int_t i=6;i<15;i++){
205     h[5]->Add(h[i]);
206     h2[5]->Add(h2[i]);
207     hPsi[i]->Draw("SAME");
208     hPsi[i]->Scale(scalA / intA * 10);
209   }
210
211   // NUA correction check V0C
212   TCanvas *can2 = new TCanvas("cV0Ccheck","cV0Ccheck");
213   hPhiPC->Draw();
214   hPsi2[5]->Draw("SAME");
215   for(Int_t i=6;i<15;i++){
216     h[5]->Add(h[i]);
217     h2[5]->Add(h2[i]);
218     hPsi2[i]->Draw("SAME");
219     hPsi2[i]->Scale(scalC /intC * 10);
220   }
221   hPsi2[5]->Scale(scalC /intC * 10);
222
223   // Flow for V0A and V0C separately
224   TCanvas *can3 = new TCanvas("cFlowComp","cFlowComp");
225   TProfile *pp = h[5]->ProfileY();
226   pp->Draw();
227   TProfile *pp2 = h2[5]->ProfileY();
228   pp2->Draw("SAME");
229
230   printf("nev (selected within 0-80%) = %i\n",p1->GetEntries());
231
232
233   // correction for resoltion
234   Float_t scaling = sqrt(res1*res3/res2);
235   pp->Scale(1./scaling);
236   printf("resolution V0A = %f\n",scaling);
237   Float_t err1_2 = eres1*eres1/res1/res1/4 +
238                    eres2*eres2/res2/res2/4 +
239                    eres3*eres3/res3/res3/4;
240   Float_t err2_2 = err1_2;
241   err1_2 /= scaling*scaling;
242   scaling = sqrt(res2*res3/res1);
243   err2_2 /= scaling*scaling;
244   pp2->Scale(1./scaling);
245   printf("resolution V0C =%f\n",scaling);
246
247   char title[100];
248   sprintf(title,"VZERO EP;p_{t} (GeV/c);v_{%i}",arm);
249   pp->SetTitle(title);
250
251   // Average V0A-V0C
252   TH1D *pAll = pp->ProjectionX();
253
254   for(Int_t i=1;i <= pAll->GetNbinsX();i++){
255        Float_t e1 = err1_2*pp->GetBinContent(i)*pp->GetBinContent(i) + pp->GetBinError(i)*pp->GetBinError(i);
256        Float_t e2 = err2_2*pp2->GetBinContent(i)*pp2->GetBinContent(i) + pp2->GetBinError(i)*pp2->GetBinError(i);
257        Float_t xval = 0,exval = 0;
258        if(e1 >0 && e2>0){
259          xval = (pp->GetBinContent(i)/e1 + pp2->GetBinContent(i)/e2)/(1/e1 + 1/e2);     
260          exval = 1./sqrt(1/e1 + 1/e2);
261        }
262        pAll->SetBinContent(i,xval);
263        pAll->SetBinError(i,exval);
264   }
265   // combined measurement
266   TCanvas *can4 = new TCanvas("cFlowCombined","cFlowCombined");
267   pAll->Draw();
268
269   char name[100];
270   sprintf(name,"out%i-%i_%i_%4.2f_%i%sv%i.root",iMin[0],iMax[0],step,pTh,charge,iTechName[itech],arm);
271   TFile *fout = new TFile(name,"RECREATE");
272
273   pAll->SetName("histo");
274   pAll->Write();
275   can1->Write();
276   can2->Write();
277   can3->Write();
278   fout->Close();
279 }