]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FLOW/macros/extractFlowVZERO.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / extractFlowVZERO.C
CommitLineData
0f25ad32 1TFile *fo = NULL;
2TF1 *fPsi = NULL;
3
4Bool_t kLib = kFALSE;
5Bool_t kVZEROrescorr = kTRUE; // VZERO resolution corrections
6Bool_t kNUAcorr = kTRUE; // NUA corrections
7Bool_t kNUOcorr = kFALSE; // Not Uniform Occupancy (TOF) corrections
8Bool_t kPIDcorr = kFALSE; // PID corrections
ad14379c 9Bool_t kFEEDcorr = kFALSE; // antiproton feed-down corrections
0f25ad32 10Bool_t kCleanMemory = kTRUE; // if kFALSE take care of the large numbers of TCanvases
11
12TH1D *hNUO[8]; // NUO corrections
13
14void LoadLib();
ad14379c 15void extractFlowVZERO(Int_t icentr,const char *type="",Int_t arm=2,Float_t pTh=0.8,Bool_t isMC=kFALSE,Int_t addbin=0);
0f25ad32 16TProfile *extractFlowVZEROsingle(Int_t icentr,Int_t spec,Int_t arm=2,Bool_t isMC=kFALSE,Float_t pTh=0.9,Int_t addbin=0,const char *nameSp="",Float_t detMin=0,Float_t detMax=1,Int_t chMin=-1,Int_t chMax=1);
17void compareTPCTOF(Int_t icentr,Int_t spec,Int_t arm=2,Float_t pTh=0.9,Int_t addbin=0);
18void plotQApid(Int_t ic,Float_t pt,Int_t addbin=0);
19
20void LoadLib(){
21 if(! kLib){
4070f709 22 gSystem->Load("libVMC");
23 gSystem->Load("libPhysics");
24 gSystem->Load("libTree");
25 gSystem->Load("libSTEERBase");
26 gSystem->Load("libANALYSIS");
27 gSystem->Load("libAOD");
28 gSystem->Load("libESD");
29 gSystem->Load("libANALYSIS");
30 gSystem->Load("libANALYSISalice");
31 gSystem->Load("libCORRFW");
32 gSystem->Load("libNetx");
33 gSystem->Load("libPWGflowBase");
34 gSystem->Load("libPWGflowTasks");
0f25ad32 35
36 kLib = kTRUE;
37 }
38}
39
ad14379c 40void extractFlowVZERO(Int_t icentr,const char *type,Int_t arm,Float_t pTh,Bool_t isMC,Int_t addbin){
0f25ad32 41 LoadLib();
42
ad14379c 43 char name[100];
44 snprintf(name,100,"AnalysisResults%s.root",type);
45 if(!fo) fo = new TFile(name);
46
0f25ad32 47 new TCanvas();
48
49 Int_t cMin[] = {0,5,10,20,30,40,50,60,70};
50 Int_t cMax[] = {5,10,20,30,40,50,60,70,80};
51
52 if(kNUOcorr){ // Compute correction for NUO in TOF
53 compareTPCTOF(icentr,0,arm,pTh,addbin);
54// compareTPCTOF(icentr,1,arm,pTh,addbin);
55// compareTPCTOF(icentr,2,arm,pTh,addbin);
56// compareTPCTOF(icentr,3,arm,pTh,addbin);
57// compareTPCTOF(icentr,4,arm,pTh,addbin);
58 }
59
60 TProfile *pAll;
61 pAll=extractFlowVZEROsingle(icentr,0,arm,0,pTh,addbin,"all",0,1);
62 pAll->SetMarkerStyle(24);
63 TProfile *pPiTOF,*pPiTPC,*pPiTPC2;
64 pPiTOF=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTOF",1,1);
65 pPiTPC=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC",0,0);
66 pPiTPC2=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC2",2,2);
67 pPiTPC->Add(pPiTPC2);
68
69 TH1D *hPi = pPiTOF->ProjectionX("hPi");
70 hPi->SetLineColor(4);
71 hPi->SetMarkerColor(4);
72 hPi->SetMarkerStyle(20);
73 for(Int_t i=1;i <=hPi->GetNbinsX();i++){
74 Float_t x = hPi->GetBinCenter(i);
75 if(x < 0.2){
76 hPi->SetBinContent(i,0);
77 hPi->SetBinError(i,0);
78 }
79 else if(x < 0.5){
80 hPi->SetBinContent(i,pPiTPC->GetBinContent(i));
81 hPi->SetBinError(i,pPiTPC->GetBinError(i));
82 }
83 else{
84 if(kNUOcorr){
85 hPi->SetBinContent(i,pPiTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
86 hPi->SetBinError(i,pPiTOF->GetBinError(i));
87 }
88 else{
89 hPi->SetBinContent(i,pPiTOF->GetBinContent(i));
90 hPi->SetBinError(i,pPiTOF->GetBinError(i));
91 }
92 }
93 }
94
95 TProfile *pElTOF,*pElTPC,*pElTPC2;
96 pElTOF=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTOF",1,1);
97 pElTPC=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC",0,0);
98 pElTPC2=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC2",2,2);
99 pElTPC->Add(pElTPC2);
100
101 TH1D *hEl = pElTOF->ProjectionX("hEl");
102 hEl->SetLineColor(6);
103 hEl->SetMarkerColor(6);
104 hEl->SetMarkerStyle(20);
105 for(Int_t i=1;i <=hEl->GetNbinsX();i++){
106 Float_t x = hEl->GetBinCenter(i);
107 if(x < 0.2){
108 hEl->SetBinContent(i,0);
109 hEl->SetBinError(i,0);
110 }
111 else if(x < 0.3){
112 hEl->SetBinContent(i,pElTPC->GetBinContent(i));
113 hEl->SetBinError(i,pElTPC->GetBinError(i));
114 }
115 else{
116 if(kNUOcorr){
117 hEl->SetBinContent(i,pElTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
118 hEl->SetBinError(i,pElTOF->GetBinError(i));
119 }
120 else{
121 hEl->SetBinContent(i,pElTOF->GetBinContent(i));
122 hEl->SetBinError(i,pElTOF->GetBinError(i));
123 }
124 }
125 }
126
127 TProfile *pKTOF,*pKTPC,*pKTPC2;
128 pKTOF=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTOF",1,1);
129 pKTPC=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC",0,0);
130 pKTPC2=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC2",2,2);
131 pKTPC->Add(pKTPC2);
132
133 TH1D *hK = pKTOF->ProjectionX("hKa");
134 hK->SetLineColor(1);
135 hK->SetMarkerColor(1);
136 hK->SetMarkerStyle(21);
137 for(Int_t i=1;i <=hK->GetNbinsX();i++){
138 Float_t x = hK->GetBinCenter(i);
139 if(x < 0.25){
140 hK->SetBinContent(i,0);
141 hK->SetBinError(i,0);
142 }
143 else if(x < 0.45){
144 hK->SetBinContent(i,pKTPC->GetBinContent(i));
145 hK->SetBinError(i,pKTPC->GetBinError(i));
146 }
147 else{
148 if(kNUOcorr){
149 hK->SetBinContent(i,pKTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
150 hK->SetBinError(i,pKTOF->GetBinError(i));
151 }
152 else{
153 hK->SetBinContent(i,pKTOF->GetBinContent(i));
154 hK->SetBinError(i,pKTOF->GetBinError(i));
155 }
156 }
157 }
158
159 TProfile *pPrTOF,*pPrTOF2,*pPrTPC,*pPrTPC2;
160 pPrTOF=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTOF",1,1,-1,-1);
161 pPrTOF2=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTOF2",1,1,-1,1);
162 pPrTPC=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTPC",0,0,-1,-1);
163 pPrTPC2=extractFlowVZEROsingle(icentr,3,arm,0,pTh,addbin,"prTPC2",2,2,-1,-1);
164 pPrTPC->Add(pPrTPC2);
165
166 TH1D *hPr = pPrTOF->ProjectionX("hPr");
167 hPr->SetLineColor(2);
168 hPr->SetMarkerColor(2);
169 hPr->SetMarkerStyle(22);
170 for(Int_t i=1;i <=hPr->GetNbinsX();i++){
171 Float_t x = hPr->GetBinCenter(i);
172 if(x < 0.3){
173 hPr->SetBinContent(i,0);
174 hPr->SetBinError(i,0);
175 }
176 else if(x < 1.0){
177 hPr->SetBinContent(i,pPrTPC->GetBinContent(i));
178 hPr->SetBinError(i,pPrTPC->GetBinError(i));
179 }
180 else{
181 if(x < 3){
182 if(kNUOcorr){
183 hPr->SetBinContent(i,pPrTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
184 hPr->SetBinError(i,pPrTOF->GetBinError(i));
185 }
186 else{
187 hPr->SetBinContent(i,pPrTOF->GetBinContent(i));
188 hPr->SetBinError(i,pPrTOF->GetBinError(i));
189 }
190 }
191 else{
192 if(kNUOcorr){
193 hPr->SetBinContent(i,pPrTOF2->GetBinContent(i) + hNUO[0]->GetBinContent(i));
194 hPr->SetBinError(i,pPrTOF2->GetBinError(i));
195 }
196 else{
197 hPr->SetBinContent(i,pPrTOF2->GetBinContent(i));
198 hPr->SetBinError(i,pPrTOF2->GetBinError(i));
199 }
200 }
201 }
202 }
203
204 pAll->Draw();
205 hPi->Draw("SAME");
206 hK->Draw("SAME");
207 hPr->Draw("SAME");
208
209
210 char name[100];
211
ad14379c 212
0f25ad32 213 // PID correction
214 if(kPIDcorr){
215 TFile *fPidTOF = new TFile("$ALICE_ROOT/PWGCF/FLOW/other/BayesianPIDcontTPCTOF.root");
216 TFile *fPidTPC = new TFile("$ALICE_ROOT/PWGCF/FLOW/other/BayesianPIDcontTPC.root");
217 // pi histos
218 sprintf(name,"Pi_IDas_Picentr%i_pth%4.2f",icentr,pTh);
219 TH1D *hPidPiPi=(TH1D *) fPidTOF->Get(name);
220 sprintf(name,"Pi_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
221 TH1D *hPidPiEl=(TH1D *) fPidTOF->Get(name);
222 sprintf(name,"Pi_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
223 TH1D *hPidPiKa=(TH1D *) fPidTOF->Get(name);
224 sprintf(name,"Pi_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
225 TH1D *hPidPiPr=(TH1D *) fPidTOF->Get(name);
226 TH1D *hPidAll = new TH1D(*hPidPiPi);
227 hPidAll->Add(hPidPiKa);
228 hPidAll->Add(hPidPiPr);
229 hPidAll->Add(hPidPiEl);
230 hPidPiPi->Divide(hPidAll);
231 hPidPiKa->Divide(hPidAll);
232 hPidPiPr->Divide(hPidAll);
233 hPidPiEl->Divide(hPidAll);
234
235 sprintf(name,"Pi_IDas_Picentr%i_pth%4.2f",icentr,pTh);
236 TH1D *hPidPiPiTPC=(TH1D *) fPidTPC->Get(name);
237 sprintf(name,"Pi_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
238 TH1D *hPidPiElTPC=(TH1D *) fPidTPC->Get(name);
239 sprintf(name,"Pi_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
240 TH1D *hPidPiKaTPC=(TH1D *) fPidTPC->Get(name);
241 sprintf(name,"Pi_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
242 TH1D *hPidPiPrTPC=(TH1D *) fPidTPC->Get(name);
243 hPidAll->Reset();
244 hPidAll->Add(hPidPiPiTPC);
245 hPidAll->Add(hPidPiKaTPC);
246 hPidAll->Add(hPidPiPrTPC);
247 hPidAll->Add(hPidPiElTPC);
248 hPidPiPiTPC->Divide(hPidAll);
249 hPidPiKaTPC->Divide(hPidAll);
250 hPidPiPrTPC->Divide(hPidAll);
251 hPidPiElTPC->Divide(hPidAll);
252
253 // K histos
254 sprintf(name,"Ka_IDas_Picentr%i_pth%4.2f",icentr,pTh);
255 TH1D *hPidKaPi=(TH1D *) fPidTOF->Get(name);
256 sprintf(name,"Ka_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
257 TH1D *hPidKaEl=(TH1D *) fPidTOF->Get(name);
258 sprintf(name,"Ka_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
259 TH1D *hPidKaKa=(TH1D *) fPidTOF->Get(name);
260 sprintf(name,"Ka_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
261 TH1D *hPidKaPr=(TH1D *) fPidTOF->Get(name);
262 hPidAll->Reset();
263 hPidAll->Add(hPidKaPi);
264 hPidAll->Add(hPidKaKa);
265 hPidAll->Add(hPidKaPr);
266 hPidAll->Add(hPidKaEl);
267 hPidKaPi->Divide(hPidAll);
268 hPidKaKa->Divide(hPidAll);
269 hPidKaPr->Divide(hPidAll);
270 hPidKaEl->Divide(hPidAll);
271
272 sprintf(name,"Ka_IDas_Picentr%i_pth%4.2f",icentr,pTh);
273 TH1D *hPidKaPiTPC=(TH1D *) fPidTPC->Get(name);
274 sprintf(name,"Ka_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
275 TH1D *hPidKaElTPC=(TH1D *) fPidTPC->Get(name);
276 sprintf(name,"Ka_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
277 TH1D *hPidKaKaTPC=(TH1D *) fPidTPC->Get(name);
278 sprintf(name,"Ka_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
279 TH1D *hPidKaPrTPC=(TH1D *) fPidTPC->Get(name);
280 hPidAll->Reset();
281 hPidAll->Add(hPidKaPiTPC);
282 hPidAll->Add(hPidKaKaTPC);
283 hPidAll->Add(hPidKaPrTPC);
284 hPidAll->Add(hPidKaElTPC);
285 hPidKaPiTPC->Divide(hPidAll);
286 hPidKaKaTPC->Divide(hPidAll);
287 hPidKaPrTPC->Divide(hPidAll);
288 hPidKaElTPC->Divide(hPidAll);
289
290 // pr histos
291 sprintf(name,"Pr_IDas_Picentr%i_pth%4.2f",icentr,pTh);
292 TH1D *hPidPrPi=(TH1D *) fPidTOF->Get(name);
293 sprintf(name,"Pr_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
294 TH1D *hPidPrEl=(TH1D *) fPidTOF->Get(name);
295 sprintf(name,"Pr_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
296 TH1D *hPidPrKa=(TH1D *) fPidTOF->Get(name);
297 sprintf(name,"Pr_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
298 TH1D *hPidPrPr=(TH1D *) fPidTOF->Get(name);
299 hPidAll->Reset();
300 hPidAll->Add(hPidPrPi);
301 hPidAll->Add(hPidPrKa);
302 hPidAll->Add(hPidPrPr);
303 hPidAll->Add(hPidPrEl);
304 hPidPrPi->Divide(hPidAll);
305 hPidPrKa->Divide(hPidAll);
306 hPidPrPr->Divide(hPidAll);
307 hPidPrEl->Divide(hPidAll);
308
309 sprintf(name,"Pr_IDas_Picentr%i_pth%4.2f",icentr,pTh);
310 TH1D *hPidPrPiTPC=(TH1D *) fPidTPC->Get(name);
311 sprintf(name,"Pr_IDas_Elcentr%i_pth%4.2f",icentr,pTh);
312 TH1D *hPidPrElTPC=(TH1D *) fPidTPC->Get(name);
313 sprintf(name,"Pr_IDas_Kacentr%i_pth%4.2f",icentr,pTh);
314 TH1D *hPidPrKaTPC=(TH1D *) fPidTPC->Get(name);
315 sprintf(name,"Pr_IDas_Prcentr%i_pth%4.2f",icentr,pTh);
316 TH1D *hPidPrPrTPC=(TH1D *) fPidTPC->Get(name);
317 hPidAll->Reset();
318 hPidAll->Add(hPidPrPiTPC);
319 hPidAll->Add(hPidPrKaTPC);
320 hPidAll->Add(hPidPrPrTPC);
321 hPidAll->Add(hPidPrElTPC);
322 hPidPrPiTPC->Divide(hPidAll);
323 hPidPrKaTPC->Divide(hPidAll);
324 hPidPrPrTPC->Divide(hPidAll);
325 hPidPrElTPC->Divide(hPidAll);
326
327 for(Int_t k=1;k <=hPi->GetNbinsX();k++){
328 Float_t pt = hPi->GetBinCenter(k);
329 Float_t xPi=hPi->GetBinContent(k)*hPidPiPi->Interpolate(pt) + hK->GetBinContent(k)*hPidPiKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidPiPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidPiEl->Interpolate(pt);
330 if(pt < 0.5)
331 xPi=hPi->GetBinContent(k)*hPidPiPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidPiKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidPiPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidPiElTPC->Interpolate(pt);
332 Float_t xKa=hPi->GetBinContent(k)*hPidKaPi->Interpolate(pt) + hK->GetBinContent(k)*hPidKaKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidKaPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidKaEl->Interpolate(pt);
333 if(pt < 0.45)
334 xKa=hPi->GetBinContent(k)*hPidKaPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidKaKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidKaPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidKaElTPC->Interpolate(pt);
335 Float_t xPr=hPi->GetBinContent(k)*hPidPrPi->Interpolate(pt) + hK->GetBinContent(k)*hPidPrKa->Interpolate(pt) + hPr->GetBinContent(k)*hPidPrPr->Interpolate(pt) + hEl->GetBinContent(k)*hPidPrEl->Interpolate(pt);
336 if(pt < 1)
337 xPr=hPi->GetBinContent(k)*hPidPrPiTPC->Interpolate(pt) + hK->GetBinContent(k)*hPidPrKaTPC->Interpolate(pt) + hPr->GetBinContent(k)*hPidPrPrTPC->Interpolate(pt) + hEl->GetBinContent(k)*hPidPrElTPC->Interpolate(pt);
338
339 hPi->SetBinContent(k,hPi->GetBinContent(k)*2 - xPi);
340 hK->SetBinContent(k,hK->GetBinContent(k)*2 - xKa);
341 hPr->SetBinContent(k,hPr->GetBinContent(k)*2 - xPr);
342 }
343 }
344
ad14379c 345
346 // antiproton Feed down
347 TProfile *pFromLambda = extractFlowVZEROsingle(icentr,11,arm,0,pTh,addbin,"pFromLambda",1,1);
348 TProfile *piFromK = extractFlowVZEROsingle(icentr,12,arm,0,pTh,addbin,"piFromK",1,1,1,1);
349 TProfile *pFromLambda2 = extractFlowVZEROsingle(icentr,11,arm,0,0.6,addbin,"pFromLambdanoPID",0,1);
350 TProfile *piFromK2 = extractFlowVZEROsingle(icentr,12,arm,0,0.6,addbin,"piFromKnoPID",0,1);
351 TProfile *piFromK3 = extractFlowVZEROsingle(icentr,12,arm,0,0.6,addbin,"piFromKnoPIDtof",1,1);
352
353 TH1D *hFeedSyst = NULL;
354
355 if(kFEEDcorr){
356 hFeedSyst = new TH1D(*hPr);
357 hFeedSyst->SetName("hFeedSyst");
358 hFeedSyst->Reset();
359 for(Int_t k=1;k <=hPr->GetNbinsX();k++){
360 Float_t contam = 3.23174e-01 * TMath::Exp(- 9.46743e-01 * hPr->GetBinCenter(k));
361 Float_t corr = contam * pFromLambda->GetBinContent(k)/(1-contam);
362 Float_t corrErr = contam * pFromLambda->GetBinError(k)/(1-contam);
363 Float_t value = hPr->GetBinContent(k)/(1-contam) - corr;
364 Float_t valueErr = hPr->GetBinError(k)/(1-contam);
365
366 hFeedSyst->SetBinContent(k,hPr->GetBinContent(k) - value);
367 hFeedSyst->SetBinContent(k,sqrt(corrErr*corrErr + valueErr*valueErr - hPr->GetBinError(k)*hPr->GetBinError(k)));
368
369 hPr->SetBinContent(k,value);
370 hPr->SetBinError(k,sqrt(corrErr*corrErr + valueErr*valueErr));
371
372 }
373 hFeedSyst->Divide(hPr);
374 }
375
0f25ad32 376 // write output
ad14379c 377 snprintf(name,100,"results%03i-%03iv%i_pTh%3.1f%s.root",cMin[icentr],cMax[icentr+addbin],arm,pTh,type);
0f25ad32 378 TFile *fout = new TFile(name,"RECREATE");
379 pAll->ProjectionX()->Write();
380 hPi->Write();
381 hK->Write();
382 hPr->Write();
383 if(isMC){
384 TH1D *pTmp = extractFlowVZEROsingle(icentr,0,arm,1,pTh,addbin,"allMC",1,1,-1,1)->ProjectionX();
385 pTmp->SetLineColor(6);
386 pTmp->SetMarkerColor(6);
387 pTmp->SetMarkerStyle(24);
388 pTmp->Write();
389 pTmp = extractFlowVZEROsingle(icentr,1,arm,1,pTh,addbin,"piMC",1,1,-1,1)->ProjectionX();
390 pTmp->SetLineColor(4);
391 pTmp->SetMarkerColor(4);
392 pTmp->SetMarkerStyle(24);
393 pTmp->Write();
394 pTmp = extractFlowVZEROsingle(icentr,2,arm,1,pTh,addbin,"kaMC",1,1,-1,1)->ProjectionX();
395 pTmp->SetLineColor(1);
396 pTmp->SetMarkerColor(1);
397 pTmp->SetMarkerStyle(25);
398 pTmp->Write();
399 pTmp = extractFlowVZEROsingle(icentr,3,arm,1,pTh,addbin,"prMC",1,1,-1,-1)->ProjectionX();
400 pTmp->SetLineColor(2);
401 pTmp->SetMarkerColor(2);
402 pTmp->SetMarkerStyle(26);
403 pTmp->Write();
404 }
ad14379c 405 extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kProf")->Write();
406 extractFlowVZEROsingle(icentr,9,arm,0,pTh,addbin,"ks",0,1,1,1)->Write();
407 extractFlowVZEROsingle(icentr,9,arm,0,pTh,addbin,"ksMy",0,1,-1,-1)->Write();
408 extractFlowVZEROsingle(icentr,10,arm,0,pTh,addbin,"lambda")->Write();
409 pFromLambda->Write();
410 piFromK->Write();
411 pFromLambda2->Write();
412 piFromK2->Write();
413 piFromK3->Write();
414 if(hFeedSyst) hFeedSyst->Write();
0f25ad32 415 fout->Close();
416}
417
418TProfile *extractFlowVZEROsingle(Int_t icentr,Int_t spec,Int_t arm,Bool_t isMC,Float_t pTh,Int_t addbin,const char *nameSp,Float_t detMin,Float_t detMax,Int_t chMin,Int_t chMax){
419 LoadLib();
420
421 pTh += 0.00001;
422
587d006a 423 // NUA correction currently are missing
424 char name[100];
0f25ad32 425 char stringa[200];
426
427 snprintf(name,100,"AnalysisResults.root");
428 if(!fo) fo = new TFile(name);
d6a1c304 429 snprintf(name,100,"contVZEROv%i",arm);
587d006a 430 TList *cont = (TList *) fo->Get(name);
afa8df58 431
0f25ad32 432 cont->ls();
433
434 Float_t xMin[5] = {icentr/*centrality bin*/,chMin/*charge*/,pTh/*prob*/,-TMath::Pi()/arm/*Psi*/,detMin/*PID mask*/};
ad14379c 435 Float_t xMax[5] = {icentr+addbin,chMax,1.0,TMath::Pi()/arm,detMax};
afa8df58 436
587d006a 437 cont->ls();
afa8df58 438
587d006a 439 TProfile *p1 = cont->At(2);
440 TProfile *p2 = cont->At(3);
441 TProfile *p3 = cont->At(4);
afa8df58 442
0f25ad32 443 TH2F *hPsi2DA = cont->At(5);
444 TH2F *hPsi2DC = cont->At(6);
445
446 TH1D *hPsiA = hPsi2DA->ProjectionY("PsiA",icentr+1,icentr+addbin+1);
447 TH1D *hPsiC = hPsi2DC->ProjectionY("PsiC",icentr+1,icentr+addbin+1);
448
449 if(!fPsi) fPsi = new TF1("fPsi","pol0",-TMath::Pi()/arm,TMath::Pi()/arm);
450 hPsiA->Fit(fPsi,"0");
451 Float_t offsetA = fPsi->GetParameter(0);
452 hPsiC->Fit(fPsi,"0");
453 Float_t offsetC = fPsi->GetParameter(0);
454
455 Int_t nbinPsi = hPsiA->GetNbinsX();
456
457 Float_t *NUAcorrA = new Float_t[nbinPsi];
458 Float_t *NUAcorrC = new Float_t[nbinPsi];
459
460 for(Int_t i=0;i < nbinPsi;i++){
461 NUAcorrA[i] = offsetA/(hPsiA->GetBinContent(i+1));
462 NUAcorrC[i] = offsetC/(hPsiC->GetBinContent(i+1));
463 }
464
afa8df58 465 Float_t res1=0,res2=0,res3=0;
466 Float_t eres1=0,eres2=0,eres3=0;
467
0f25ad32 468 for(Int_t i = icentr; i <= icentr+addbin;i++){
469 if(p1->GetBinError(i+1)){
470 eres1 += 1./p1->GetBinError(i+1)/p1->GetBinError(i+1);
471 res1 += p1->GetBinContent(i+1)/p1->GetBinError(i+1)/p1->GetBinError(i+1);
472 }
473 if(p2->GetBinError(i+1)){
474 eres2 += 1./p2->GetBinError(i+1)/p2->GetBinError(i+1);
475 res2 += p2->GetBinContent(i+1)/p2->GetBinError(i+1)/p2->GetBinError(i+1);
476 }
477 if(p3->GetBinError(i+1)){
478 eres3 += 1./p3->GetBinError(i+1)/p3->GetBinError(i+1);
479 res3 += p3->GetBinContent(i+1)/p3->GetBinError(i+1)/p3->GetBinError(i+1);
afa8df58 480 }
afa8df58 481 }
0f25ad32 482
587d006a 483 res1 /= eres1;
484 res2 /= eres2;
485 res3 /= eres3;
486
0f25ad32 487 eres1 = sqrt(1./eres1);
488 eres2 = sqrt(1./eres2);
489 eres3 = sqrt(1./eres3);
490
587d006a 491 AliFlowVZEROResults *a = (AliFlowVZEROResults *) cont->At(0);
492 AliFlowVZEROResults *b = (AliFlowVZEROResults *) cont->At(1);
0f25ad32 493 TProfile *pp,*pp2;
494 if(kNUAcorr){ // with NUA corrections
495 pp = a->GetV2reweight(spec,xMin,xMax,3,NUAcorrA);
496 pp2 = b->GetV2reweight(spec,xMin,xMax,3,NUAcorrC);
497 }
498 else{
499 pp = a->GetV2(spec,xMin,xMax);
500 pp2 = b->GetV2(spec,xMin,xMax);
501 }
587d006a 502
afa8df58 503 Float_t scaling = sqrt(res1*res3/res2);
0f25ad32 504 if(kVZEROrescorr){
505 pp->Scale(1./scaling);
506 }
507
afa8df58 508 Float_t err1_2 = eres1*eres1/res1/res1/4 +
587d006a 509 eres2*eres2/res2/res2/4 +
510 eres3*eres3/res3/res3/4;
afa8df58 511 Float_t err2_2 = err1_2;
512 err1_2 /= scaling*scaling;
0f25ad32 513 printf("resolution V0A = %f +/- %f\n",scaling,err1_2);
afa8df58 514 scaling = sqrt(res2*res3/res1);
515 err2_2 /= scaling*scaling;
0f25ad32 516 if(kVZEROrescorr){
517 pp2->Scale(1./scaling);
518 }
519 printf("resolution V0C =%f +/- %f\n",scaling,err2_2);
afa8df58 520
587d006a 521 pp->SetName("V0A");
522 pp2->SetName("V0C");
afa8df58 523
0f25ad32 524 if(! kCleanMemory){
525 new TCanvas();
526 pp->Draw();
527 pp2->Draw("SAME");
528 }
529
530 TProfile *pData = new TProfile(*pp);
531 pData->Add(pp2);
532 snprintf(stringa,100,"%sData",nameSp);
533 pData->SetName(stringa);
afa8df58 534
0f25ad32 535 TProfile *pMc = NULL;
536
537 TProfile *ppMC;
538 TProfile *ppMC2;
539
540 if(isMC){
d6a1c304 541 snprintf(name,100,"contVZEROmc");
587d006a 542 cont = (TList *) fo->Get(name);
0f25ad32 543 cont->ls();
544 if(arm == 2){
545 AliFlowVZEROResults *c = (AliFlowVZEROResults *) cont->At(0);
546 if(! kCleanMemory) c->GetV2(spec,xMin,xMax)->Draw("SAME");
547 }
548 AliFlowVZEROResults *cA;
549 if(fo->Get("contVZEROv2")) cA = (AliFlowVZEROResults *) cont->At(1+2*(arm==3));
550 else cA = (AliFlowVZEROResults *) cont->At(0);
551 AliFlowVZEROResults *cC;
552 if(fo->Get("contVZEROv2")) cC = (AliFlowVZEROResults *) cont->At(2+2*(arm==3));
553 else cC = (AliFlowVZEROResults *) cont->At(1);
554
555 TProfile *p1mc = cont->At(5+3*(arm==3));
556 TProfile *p2mc = cont->At(6+3*(arm==3));
557 TProfile *p3mc = cont->At(7+3*(arm==3));
558
559 Float_t resMC1=0,resMC2=0,resMC3=0;
560 Float_t eresMC1=0,eresMC2=0,eresMC3=0;
561
562 for(Int_t i = icentr; i <= icentr+addbin;i++){
563 if(p1mc->GetBinError(i+1)){
564 eresMC1 += 1./p1mc->GetBinError(i+1)/p1mc->GetBinError(i+1);
565 resMC1 += p1mc->GetBinContent(i+1)/p1mc->GetBinError(i+1)/p1mc->GetBinError(i+1);
566 }
567 if(p2mc->GetBinError(i+1)){
568 eresMC2 += 1./p2mc->GetBinError(i+1)/p2mc->GetBinError(i+1);
569 resMC2 += p2mc->GetBinContent(i+1)/p2mc->GetBinError(i+1)/p2mc->GetBinError(i+1);
570 }
571 if(p3mc->GetBinError(i+1)){
572 eresMC3 += 1./p3mc->GetBinError(i+1)/p3mc->GetBinError(i+1);
573 resMC3 += p3mc->GetBinContent(i+1)/p3mc->GetBinError(i+1)/p3mc->GetBinError(i+1);
574 }
575 }
576
577 resMC1 /= eresMC1;
578 resMC2 /= eresMC2;
579 resMC3 /= eresMC3;
580
581 eresMC1 = sqrt(1./eresMC1);
582 eresMC2 = sqrt(1./eresMC2);
583 eresMC3 = sqrt(1./eresMC3);
584
585 ppMC = cA->GetV2(spec,xMin,xMax);
586 ppMC2 = cC->GetV2(spec,xMin,xMax);
587
588 scaling = sqrt(resMC1*resMC3/resMC2);
589 ppMC->Scale(1./scaling);
590
591 err1_2 = eresMC1*eresMC1/resMC1/resMC1/4 +
592 eresMC2*eresMC2/resMC2/resMC2/4 +
593 eresMC3*eresMC3/resMC3/resMC3/4;
594 err2_2 = err1_2;
595 err1_2 /= scaling*scaling;
596 printf("resolution V0A (MC) = %f +/- %f\n",scaling,err1_2);
597 scaling = sqrt(resMC2*resMC3/resMC1);
598 err2_2 /= scaling*scaling;
599 ppMC2->Scale(1./scaling);
600 printf("resolution V0C (MC) =%f +/- %f\n",scaling,err2_2);
601
602 ppMC->SetName("V0Amc");
603 ppMC2->SetName("V0Cmc");
604
605 if(! kCleanMemory){
606 ppMC->Draw("SAME");
607 ppMC2->Draw("SAME");
608 }
609
610 pMc = new TProfile(*ppMC);
611 pMc->Add(ppMC2);
612 snprintf(stringa,100,"%sMC",nameSp);
613 pMc->SetName(stringa);
614 pMc->SetLineColor(2);
afa8df58 615 }
0f25ad32 616
617 if(! kCleanMemory){
618 new TCanvas();
619 pData->Draw();
620 }
621 if(pMc && !kCleanMemory){
622 pMc->Draw("SAME");
623 TH1D *hData = pData->ProjectionX();
624 TH1D *hMc = pMc->ProjectionX();
625 hData->Divide(hMc);
626 new TCanvas();
627 hData->Draw();
628 }
629
630 delete[] NUAcorrA;
631 delete[] NUAcorrC;
632
633 if(kCleanMemory){
634 if(pp) delete pp;
635 if(pp2) delete pp2;
636 if(ppMC) delete ppMC;
637 if(ppMC2) delete ppMC2;
638 }
639
640 delete cont;
641
642 if(isMC) return pMc;
643 return pData;
644}
645void compareTPCTOF(Int_t icentr,Int_t spec,Int_t arm,Float_t pTh,Int_t addbin){
646 LoadLib();
647 Float_t ptMaxFit[8] = {20,0.7,0.55,1.2,2,2,2,2};
648
649 TProfile *pTPC = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TPC",0,0,-1,-1+2*(spec!=3)); // TPC && !TOF (TPC PID)
650 TProfile *pTPC2 = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TPCextra",2,2,-1,-1+2*(spec!=3)); //TPC && TOF (TPC PID)
651
652 TProfile *pTOF = extractFlowVZEROsingle(icentr,spec,arm,0,pTh,addbin,"TOF",1,1,-1,-1+2*(spec!=3)); //TPC && TOF (TPC&TOF PID)
653
654 if(spec > 0) pTPC->Add(pTPC2);
655 else pTPC->Add(pTOF);
656
657 char name[100];
658 snprintf(name,100,"NUO_%i",spec);
659
660 hNUO[spec] = pTPC->ProjectionX(name);
661 hNUO[spec]->Add(pTOF,-1);
662
663 if(spec && hNUO[0]) hNUO[spec]->Add(hNUO[0],-1);
664
665 if(! kCleanMemory){
666 new TCanvas();
667 pTPC->Draw();
668 pTOF->Draw("SAME");
669
670 new TCanvas();
671 pTPC->SetLineColor(1);
672 pTOF->SetLineColor(2);
673 }
674
675 hNUO[spec]->Draw();
676
677 if(kCleanMemory){
678 if(pTPC) delete pTPC;
679 if(pTPC2) delete pTPC2;
680 if(pTOF) delete pTOF;
681 }
682}
683
684void plotQApid(Int_t ic,Float_t pt,Int_t addbin){
685 char name[100];
686 char stringa[200];
687 LoadLib();
688
689 snprintf(name,100,"AnalysisResults.root");
690 if(!fo) fo = new TFile(name);
691 snprintf(name,100,"contVZEROv2");
692 TList *cont = (TList *) fo->Get(name);
693 AliFlowVZEROResults *pidqa = cont->FindObject("qaPID");
694 Float_t xval[2] = {2.,pt+0.00001};
695 Float_t xval2[2] = {2.+addbin,pt+0.00001};
696
697 TProfile *proTPCpi = pidqa->GetV2(0,xval,xval2);
698 TProfile *proTOFpi = pidqa->GetV2(1,xval,xval2);
699 TProfile *proTPCka = pidqa->GetV2(2,xval,xval2);
700 TProfile *proTOFka = pidqa->GetV2(3,xval,xval2);
701 TProfile *proTPCpr = pidqa->GetV2(4,xval,xval2);
702 TProfile *proTOFpr = pidqa->GetV2(5,xval,xval2);
703
704 proTPCpi->SetName("hPiTPC");
705 proTOFpi->SetName("hPiTOF");
706 proTPCka->SetName("hKaTPC");
707 proTOFka->SetName("hKaTOF");
708 proTPCpr->SetName("hPrTPC");
709 proTOFpr->SetName("hPrTOF");
710
711 proTPCpi->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{#pi} (a.u)");
712 proTOFpi->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{#pi} (ps)");
713 proTPCka->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{K} (a.u)");
714 proTOFka->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{K} (ps)");
715 proTPCpr->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{p} (a.u)");
716 proTOFpr->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{p} (ps)");
717
718 TCanvas *c = new TCanvas("cVproj","cVproj");
719 c->Divide(2,3);
720 c->cd(1);
721 proTPCpi->Draw();
722 c->cd(2);
723 proTOFpi->Draw();
724 c->cd(3);
725 proTPCka->Draw();
726 c->cd(4);
727 proTOFka->Draw();
728 c->cd(5);
729 proTPCpr->Draw();
730 c->cd(6);
731 proTOFpr->Draw();
afa8df58 732}