]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/macros/readBalanceFunction.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update of macro used for extracting Balance Function widths from data
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / readBalanceFunction.C
CommitLineData
b7ce3d3f 1const TString gBFAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
2
3const Int_t nrOfCentralities = 9;
4const Double_t centralityArray[nrOfCentralities+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.}; // in centrality percentile (0-5,5-10,10-20,20-30,...,70-80)
5//const Double_t centralityArray[nrOfCentralities+1] = {0.,1.,2.,3.,4.,6.,10.,20.,30.,80.}; // in centrality percentile (0-5,5-10,10-20,20-30,...,70-80)
6const Double_t cent[nrOfCentralities] = {382.8,329.7,260.5,186.4,128.9,85.,52.8,30.,15.8}; // hard coded at the moment for centrality percentiles
7const Double_t centE[nrOfCentralities] = {3.1,4.6,4.4,3.9,3.3,2.6,2.0,1.3,0.6}; // (0-5,5-10,10-20,20-30,...,70-80)
8
ab34c777 9void readBalanceFunction(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root",Int_t fStartBinBFWidth = 3, Int_t fRebin = 2,Int_t fStartBinBFWidthPhi = 2, Int_t fRebinPhi = 2,TString centEst = "V0M",Double_t etaWindow = -1) {
b7ce3d3f 10 // Macro to read the output of the BF analysis: MW: CHANGE THIS!!!!
11 //i) Prints and draws the final BF output
12 //ii) Plots the QA part of the analysis
13 //iii) store BF in output file
14 //Author: Panos.Christakoglou@cern.ch, m.weber@cern.ch
15 //Loading the needed libraries
16 gSystem->Load("libProofPlayer.so");
17 gSystem->Load("libANALYSIS.so");
18 gSystem->Load("libANALYSISalice.so");
19 gSystem->Load("libPWGCFebye.so");
20
21 //Draw BF
ab34c777 22 drawBF(bHistos,inFile, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,centEst, "", etaWindow);
b7ce3d3f 23
24
25 //Merge the output
26 //mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");
27}
28
29//___________________________________________________________//
ab34c777 30void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", Int_t fStartBinBFWidth = 1, Int_t fRebin = 1, Int_t fStartBinBFWidthPhi = 1, Int_t fRebinPhi = 1, TString centEst = "V0M",TString extraString = "", Double_t etaWindow = -1) {
b7ce3d3f 31 //Function to draw the BF objects and write them into the output file
32
33 Int_t maximumCanvases = 13;
34 Int_t iCanvas = 0;
35 Int_t iList = -1;
36 TCanvas *cQA[13][10];
37 TCanvas *cQAV0M = new TCanvas("cQAV0M","V0M multiplicities");
38 cQAV0M->Divide(4,3);
39 TCanvas *cQARef = new TCanvas("cQARef","reference track multiplicities");
40 cQARef->Divide(4,3);
41 TCanvas *cBF[13][10];
42 TCanvas *cBFS[13][10];
43
44 // get the file
45 TFile *f = TFile::Open(inFile.Data());
46 if(!f) {
47 Printf("File not found!!!");
48 break;
49 }
50
51 // get the BF output directory
52 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis"));
53 if(!dir) {
54 Printf("Output directory not found!!!");
55 break;
56 }
57
58 // loop over all lists and plot the BF and QA
59 TList *list = NULL;
60 TString listName;
61 TIter nextkey( dir->GetListOfKeys() );
62 TKey *key;
63
64 AliBalance *bf[13][7];
65 AliBalance *bfs[13][7];
66 TH1D *gbf[13][10][7];
67 TH1D *gbfs[13][10][7];
68
69 for(Int_t i = 0; i < 13; i++){
70 for(Int_t j = 0; j < 10; j++){
71 for(Int_t k = 0; k < 7; k++){
72 gbf[i][j][k] = NULL;
73 gbfs[i][j][k] = NULL;
74 }
75 }
76 }
77
78 TH2D *fHistP[13][7]; //N+
79 TH2D *fHistN[13][7]; //N-
80 TH2D *fHistPN[13][7]; //N+-
81 TH2D *fHistNP[13][7]; //N-+
82 TH2D *fHistPP[13][7]; //N++
83 TH2D *fHistNN[13][7]; //N--
84
85 TH2D *fHistPS[13][7]; //N+
86 TH2D *fHistNS[13][7]; //N-
87 TH2D *fHistPNS[13][7]; //N+-
88 TH2D *fHistNPS[13][7]; //N-+
89 TH2D *fHistPPS[13][7]; //N++
90 TH2D *fHistNNS[13][7]; //N--
91
92 Double_t WM[13][10]; // weighted mean for eta (recalculated from fStartBin)
93 Double_t WME[13][10]; // error
94 Double_t WMS[13][10]; // weighted mean for eta (recalculated from fStartBin) (shuffled)
95 Double_t WMSE[13][10]; // error (shuffled)
96
ab34c777 97 Double_t WMP[13][10]; // weighted mean for phi (recalculated from fStartBinPhi)
b7ce3d3f 98 Double_t WMPE[13][10]; // error
ab34c777 99 Double_t WMPS[13][10]; // weighted mean for phi (recalculated from fStartBin) (shuffled)
b7ce3d3f 100 Double_t WMPSE[13][10]; // error (shuffled)
101
102 Double_t integ[13][10]; // integral for eta (calculated from bin 1)
103 Double_t integE[13][10]; // error
104 Double_t integS[13][10]; // integral for eta (calculated from bin 1) (shuffled)
105 Double_t integSE[13][10]; // error (shuffled)
106
107 Double_t integP[13][10]; // integral for phi (calculated from bin 1)
108 Double_t integPE[13][10]; // error
109 Double_t integPS[13][10]; // integral for phi (calculated from bin 1) (shuffled)
110 Double_t integPSE[13][10]; // error (shuffled)
111
112
113 while ( (key = (TKey*)nextkey())) {
114
115 list = (TList*)key->ReadObj();
116 listName = TString(list->GetName());
117
118 cout<<"Processing list "<<listName<<endl;
119
120
121 // ----------------------------------------------------
122 // plot QA histograms
123 if(listName.Contains("QA")){
124
125 if(iList<13) iList++;
126 else{
127 cerr<<"TOO MANY LISTS!!!"<<endl;
128 return;
129 }
130
131 cQA[iList][iCanvas] = new TCanvas(listName,listName);
132 cQA[iList][iCanvas]->Divide(4,3);
133
134 cQA[iList][iCanvas]->cd(1);
135 TH1F* histEventStats = (TH1F*)list->FindObject("fHistEventStats");
136 if(histEventStats){
137 histEventStats->SetFillColor(9);
138 histEventStats->Draw();
139 }
140
141 cQA[iList][iCanvas]->cd(2);
142 TH1F* histTriggerStats = (TH1F*)list->FindObject("fHistTriggerStats");
143 if(histTriggerStats){
144 histTriggerStats->SetFillColor(9);
145 histTriggerStats->Draw();
146 }
147
148 cQA[iList][iCanvas]->cd(3);
149 TH1F* histTrackStats = (TH1F*)list->FindObject("fHistTrackStats");
150 if(histTrackStats){
151 histTrackStats->SetFillColor(9);
152 histTrackStats->Draw();
153 }
154
155 cQA[iList][iCanvas]->cd(4);
156 TH1F* histCentStats = (TH1F*)list->FindObject("fHistCentStats");
157 if(histCentStats){
158 histCentStats->SetFillColor(9);
159 histCentStats->Draw("colz");
160 }
161
162
163
164 cQA[iList][iCanvas]->cd(5);
165 TH1F* histVx = (TH1F*)list->FindObject("fHistVx");
166 if(histVx){
167 histVx->SetFillColor(9);
168 histVx->Draw();
169 }
170 cQA[iList][iCanvas]->cd(6);
171 TH1F* histVy = (TH1F*)list->FindObject("fHistVy");
172 if(histVy){
173 histVy->SetFillColor(9);
174 histVy->Draw();
175 }
176
177 cQA[iList][iCanvas]->cd(7);
178 TH1F* histVz = (TH1F*)list->FindObject("fHistVz");
179 if(histVz){
180 histVz->SetFillColor(9);
181 histVz->Draw();
182 }
183
184 cQA[iList][iCanvas]->cd(8);
185 cQA[iList][iCanvas]->cd(8)->SetLogz();
186 TH2F* histDCA = (TH2F*)list->FindObject("fHistDCA");
187 if(histDCA) histDCA->Draw("colz");
188
189
190 cQA[iList][iCanvas]->cd(9);
191 cQA[iList][iCanvas]->cd(9)->SetLogz();
192 TH2F* histClus = (TH2F*)list->FindObject("fHistClus");
193 if(histClus) histClus->Draw("colz");
194
195 cQA[iList][iCanvas]->cd(10);
196 TH1F* histPt = (TH1F*)list->FindObject("fHistPt");
197 if(histPt){
198 histPt->SetFillColor(9);
199 histPt->Draw();
200 }
201
202 cQA[iList][iCanvas]->cd(11);
203 TH1F* histEta = (TH1F*)list->FindObject("fHistEta");
204 if(histEta){
205 histEta->SetFillColor(9);
206 histEta->Draw();
207 }
208
209 cQA[iList][iCanvas]->cd(12);
210 TH1F* histPhi = (TH1F*)list->FindObject("fHistPhi");
211 if(histPhi){
212 histPhi->SetFillColor(9);
213 histPhi->Draw();
214 }
215
216 // centrality estimator QA
217 cQAV0M->cd(iCanvas+1);
218 cQAV0M->cd(iCanvas+1)->SetLogz();
219 TH1F* histV0M = (TH1F*)list->FindObject("fHistV0M");
220 if(histV0M){
221 histV0M->Draw("colz");
222 }
223
224 cQARef->cd(iCanvas+1);
225 cQARef->cd(iCanvas+1)->SetLogz();
226 TH1F* histRef = (TH1F*)list->FindObject("fHistRefTracks");
227 if(histRef){
228 histRef->Draw("colz");
229 }
230 }
231 // ----------------------------------------------------
232
233 // ----------------------------------------------------
234 // calculate and plot BF
ab34c777 235 if(listName.Contains("BF_")&&listName.Contains(centEst.Data())){
b7ce3d3f 236
237 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
238
239 cBF[iList][iCanvas] = new TCanvas(Form("cBF_%d_%d",iList,iCanvas),Form("cBF_%d_%d",iList,iCanvas));
240 cBF[iList][iCanvas]->Divide(3,3);
241
242 }
243
244
245 for(Int_t a = 0; a < 7; a++){
246
247 cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;
248
249 // create the BF object
250 bf[iList][a] = new AliBalance();
251
252 fHistP[iList][a] = (TH2D*)list->FindObject(Form("fHistP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
253 fHistN[iList][a] = (TH2D*)list->FindObject(Form("fHistN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
254 fHistPP[iList][a] = (TH2D*)list->FindObject(Form("fHistPP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
255 fHistPN[iList][a] = (TH2D*)list->FindObject(Form("fHistPN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
256 fHistNP[iList][a] = (TH2D*)list->FindObject(Form("fHistNP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
257 fHistNN[iList][a] = (TH2D*)list->FindObject(Form("fHistNN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
258
259 // rebin histograms (be careful with divider!)
260 if(a==6){
ab34c777 261 fHistP[iList][a]->RebinY(fRebinPhi);
262 fHistN[iList][a]->RebinY(fRebinPhi);
263 fHistPP[iList][a]->RebinY(fRebinPhi);
264 fHistPN[iList][a]->RebinY(fRebinPhi);
265 fHistNP[iList][a]->RebinY(fRebinPhi);
266 fHistNN[iList][a]->RebinY(fRebinPhi);
b7ce3d3f 267 }
268 else{
269 fHistP[iList][a]->RebinY(fRebin);
270 fHistN[iList][a]->RebinY(fRebin);
271 fHistPP[iList][a]->RebinY(fRebin);
272 fHistPN[iList][a]->RebinY(fRebin);
273 fHistNP[iList][a]->RebinY(fRebin);
274 fHistNN[iList][a]->RebinY(fRebin);
275 }
276
277 fHistP[iList][a]->SetName(Form("%s_%d_%d",fHistP[iList][a]->GetName(),iList,iCanvas));
278 fHistN[iList][a]->SetName(Form("%s_%d_%d",fHistN[iList][a]->GetName(),iList,iCanvas));
279 fHistPP[iList][a]->SetName(Form("%s_%d_%d",fHistPP[iList][a]->GetName(),iList,iCanvas));
280 fHistPN[iList][a]->SetName(Form("%s_%d_%d",fHistPN[iList][a]->GetName(),iList,iCanvas));
281 fHistNP[iList][a]->SetName(Form("%s_%d_%d",fHistNP[iList][a]->GetName(),iList,iCanvas));
282 fHistNN[iList][a]->SetName(Form("%s_%d_%d",fHistNN[iList][a]->GetName(),iList,iCanvas));
283
284 // set histograms in AliBalance object
285 bf[iList][a]->SetHistNp(a, fHistP[iList][a]);
286 bf[iList][a]->SetHistNn(a, fHistN[iList][a]);
287 bf[iList][a]->SetHistNpp(a, fHistPP[iList][a]);
288 bf[iList][a]->SetHistNpn(a, fHistPN[iList][a]);
289 bf[iList][a]->SetHistNnp(a, fHistNP[iList][a]);
290 bf[iList][a]->SetHistNnn(a, fHistNN[iList][a]);
291
292 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
293
294 gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow);
295 gbf[iList][iCanvas][a]->SetName(Form("BF_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
296
297 cBF[iList][iCanvas]->cd(a+1);
298 gbf[iList][iCanvas][a]->SetMarkerStyle(20);
299
300 if(!bHistos){
301 gbf[iList][iCanvas][a]->DrawCopy("AP");
302 if(a==1){
ab34c777 303 GetWeightedMean(gbf[iList][iCanvas][a],fStartBinBFWidth,WM[iList][iCanvas],WME[iList][iCanvas]); // for eta recalculate width
b7ce3d3f 304 GetIntegral(gbf[iList][iCanvas][a],integ[iList][iCanvas],integE[iList][iCanvas]);
305 }
306 else if(a==6){
ab34c777 307 GetWeightedMean(gbf[iList][iCanvas][a],fStartBinBFWidthPhi,WMP[iList][iCanvas],WMPE[iList][iCanvas]); // for phi calculate width
b7ce3d3f 308 GetIntegral(gbf[iList][iCanvas][a],integP[iList][iCanvas],integPE[iList][iCanvas]);
309 }
310 }
311 else{
312 fHistPN[iList][a]->SetLineColor(2);
313 fHistPN[iList][a]->ProjectionY(Form("pn%d",a))->DrawCopy();
314 fHistPP[iList][a]->SetLineColor(1);
315 fHistPP[iList][a]->ProjectionY(Form("pp%d",a))->DrawCopy("same");
316 fHistNP[iList][a]->SetLineColor(4);
317 fHistNP[iList][a]->ProjectionY(Form("np%d",a))->DrawCopy("same");
318 fHistNN[iList][a]->SetLineColor(8);
319 fHistNN[iList][a]->ProjectionY(Form("nn%d",a))->DrawCopy("same");
320 }
321 }
322 }
323 }
324 // ----------------------------------------------------
325
326 // ----------------------------------------------------
327 // calculate and plot BF (shuffled)
ab34c777 328 if(listName.Contains("BFShuffled")&&listName.Contains(centEst.Data())&&listName.Contains(extraString.Data())){
b7ce3d3f 329
330 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
331
332 cBFS[iList][iCanvas] = new TCanvas(Form("Shuffled_%d_%d",iList,iCanvas),Form("Shuffled_%d_%d",iList,iCanvas));
333 cBFS[iList][iCanvas]->Divide(3,3);
334
335 }
336
337 for(Int_t a = 0; a < 7; a++){
338
339 // create the BF object
340 bfs[iList][a] = new AliBalance();
341
342 fHistPS[iList][a] = (TH2D*)list->FindObject(Form("fHistP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
343 fHistNS[iList][a] = (TH2D*)list->FindObject(Form("fHistN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
344 fHistPPS[iList][a] = (TH2D*)list->FindObject(Form("fHistPP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
345 fHistPNS[iList][a] = (TH2D*)list->FindObject(Form("fHistPN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
346 fHistNPS[iList][a] = (TH2D*)list->FindObject(Form("fHistNP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
347 fHistNNS[iList][a] = (TH2D*)list->FindObject(Form("fHistNN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
348
349 // rebin histograms (be careful with divider!)
350 if(a==6){
ab34c777 351 fHistPS[iList][a]->RebinY(fRebinPhi);
352 fHistNS[iList][a]->RebinY(fRebinPhi);
353 fHistPPS[iList][a]->RebinY(fRebinPhi);
354 fHistPNS[iList][a]->RebinY(fRebinPhi);
355 fHistNPS[iList][a]->RebinY(fRebinPhi);
356 fHistNNS[iList][a]->RebinY(fRebinPhi);
b7ce3d3f 357 }
358 else{
359 fHistPS[iList][a]->RebinY(fRebin);
360 fHistNS[iList][a]->RebinY(fRebin);
361 fHistPPS[iList][a]->RebinY(fRebin);
362 fHistPNS[iList][a]->RebinY(fRebin);
363 fHistNPS[iList][a]->RebinY(fRebin);
364 fHistNNS[iList][a]->RebinY(fRebin);
365 }
366
367 fHistPS[iList][a]->SetName(Form("%s_%d_%d",fHistPS[iList][a]->GetName(),iList,iCanvas));
368 fHistNS[iList][a]->SetName(Form("%s_%d_%d",fHistNS[iList][a]->GetName(),iList,iCanvas));
369 fHistPPS[iList][a]->SetName(Form("%s_%d_%d",fHistPPS[iList][a]->GetName(),iList,iCanvas));
370 fHistPNS[iList][a]->SetName(Form("%s_%d_%d",fHistPNS[iList][a]->GetName(),iList,iCanvas));
371 fHistNPS[iList][a]->SetName(Form("%s_%d_%d",fHistNPS[iList][a]->GetName(),iList,iCanvas));
372 fHistNNS[iList][a]->SetName(Form("%s_%d_%d",fHistNNS[iList][a]->GetName(),iList,iCanvas));
373
374 // set histograms in AliBalance object
375 bfs[iList][a]->SetHistNp(a, fHistPS[iList][a]);
376 bfs[iList][a]->SetHistNn(a, fHistNS[iList][a]);
377 bfs[iList][a]->SetHistNpp(a, fHistPPS[iList][a]);
378 bfs[iList][a]->SetHistNpn(a, fHistPNS[iList][a]);
379 bfs[iList][a]->SetHistNnp(a, fHistNPS[iList][a]);
380 bfs[iList][a]->SetHistNnn(a, fHistNNS[iList][a]);
381
382 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
383
384 gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow);
385 gbfs[iList][iCanvas][a]->SetName(Form("BFS_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
386
387 cBFS[iList][iCanvas]->cd(a+1);
388 gbfs[iList][iCanvas][a]->SetMarkerStyle(20);
389 if(!bHistos){
390 gbfs[iList][iCanvas][a]->DrawCopy("AP");
391 if(a==1){
392 GetWeightedMean(gbfs[iList][iCanvas][a],fStartBinBFWidth,WMS[iList][iCanvas],WMSE[iList][iCanvas]);
393 GetIntegral(gbfs[iList][iCanvas][a],integS[iList][iCanvas],integSE[iList][iCanvas]);
394 }
395 else if(a==6){
ab34c777 396 GetWeightedMean(gbfs[iList][iCanvas][a],fStartBinBFWidthPhi,WMPS[iList][iCanvas],WMPSE[iList][iCanvas]); // for phi calculate width
b7ce3d3f 397 GetIntegral(gbfs[iList][iCanvas][a],integPS[iList][iCanvas],integPSE[iList][iCanvas]);
398 }
399 }
400 else{
401 fHistPNS[iList][a]->SetLineColor(2);
402 fHistPNS[iList][a]->ProjectionY(Form("pns%d",a))->DrawCopy();
403 fHistPPS[iList][a]->SetLineColor(1);
404 fHistPPS[iList][a]->ProjectionY(Form("pps%d",a))->DrawCopy("same");
405 fHistNPS[iList][a]->SetLineColor(4);
406 fHistNPS[iList][a]->ProjectionY(Form("nps%d",a))->DrawCopy("same");
407 fHistNNS[iList][a]->SetLineColor(8);
408 fHistNNS[iList][a]->ProjectionY(Form("nns%d",a))->DrawCopy("same");
409 }
410 }
411 }
412 }
413 // ----------------------------------------------------
414 }
415
416 // for BF calculation create also graphs with weighted mean for eta
417 if(iCanvas == 0) return;
418
419 TGraphErrors *gWM[13];
420 TGraphErrors *gWMS[13];
421 TGraphErrors *gWMP[13];
422 TGraphErrors *gWMPS[13];
423 TGraphErrors *ginteg[13];
424 TGraphErrors *gintegS[13];
425 TGraphErrors *gintegP[13];
426 TGraphErrors *gintegPS[13];
427 for(Int_t i = 0; i < 13; i++){
428 gWM[i] = NULL;
429 gWMS[i] = NULL;
430 ginteg[i] = NULL;
431 gintegS[i] = NULL;
432 gWMP[i] = NULL;
433 gWMPS[i] = NULL;
434 gintegP[i] = NULL;
435 gintegPS[i] = NULL;
436 }
437
438
439 if(!bHistos){
440 for(Int_t i = 0; i < iList+1; i++){
441 gWM[i] = new TGraphErrors(iCanvas,cent,WM[i],centE,WME[i]);
442 if(iList==1) gWM[i]->SetName("gCentrality");
443 else gWM[i]->SetName(Form("gCentrality_%d",i));
444 gWMS[i] = new TGraphErrors(iCanvas,cent,WMS[i],centE,WMSE[i]);
445 if(iList==1) gWMS[i]->SetName("gCentralityS");
446 else gWMS[i]->SetName(Form("gCentralityS_%d",i));
447
448 gWMP[i] = new TGraphErrors(iCanvas,cent,WMP[i],centE,WMPE[i]);
449 if(iList==1) gWMP[i]->SetName("gCentralityPhi");
450 else gWMP[i]->SetName(Form("gCentralityPhi_%d",i));
451 gWMPS[i] = new TGraphErrors(iCanvas,cent,WMPS[i],centE,WMPSE[i]);
452 if(iList==1) gWMPS[i]->SetName("gCentralityPhiS");
453 else gWMPS[i]->SetName(Form("gCentralityPhiS_%d",i));
454
455 ginteg[i] = new TGraphErrors(iCanvas,cent,integ[i],centE,integE[i]);
456 if(iList==1) ginteg[i]->SetName("gIntegral");
457 else ginteg[i]->SetName(Form("gIntegral_%d",i));
458 gintegS[i] = new TGraphErrors(iCanvas,cent,integS[i],centE,integSE[i]);
459 if(iList==1) gintegS[i]->SetName("gIntegralS");
460 else gintegS[i]->SetName(Form("gIntegralS_%d",i));
461
462 gintegP[i] = new TGraphErrors(iCanvas,cent,integP[i],centE,integPE[i]);
463 if(iList==1) gintegP[i]->SetName("gIntegraPhil");
464 else gintegP[i]->SetName(Form("gIntegralPhi_%d",i));
465 gintegPS[i] = new TGraphErrors(iCanvas,cent,integPS[i],centE,integPSE[i]);
466 if(iList==1) gintegPS[i]->SetName("gIntegralPhiS");
467 else gintegPS[i]->SetName(Form("gIntegralPhiS_%d",i));
468 }
469 }
470
ab34c777 471 TFile *fOut = TFile::Open(Form("Histograms_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
b7ce3d3f 472 fOut->cd();
473 for(Int_t i = 0; i < iList+1; i++){
ab34c777 474 cout<<"PROCESS LIST "<<i<<" NOW!"<<endl;
b7ce3d3f 475 for(Int_t a = 0; a < 7; a++){
ab34c777 476 cout<<"PROCESS VARIABLE "<<a<<" NOW!"<<endl;
b7ce3d3f 477
478 if(fHistPN[i][a]){
479 (fHistPN[i][a]->ProjectionY(Form("hPN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
480 (fHistPP[i][a]->ProjectionY(Form("hPP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
481 (fHistNP[i][a]->ProjectionY(Form("hNP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
482 (fHistNN[i][a]->ProjectionY(Form("hNN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
483 (fHistP[i][a]->ProjectionY(Form("hP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
484 (fHistN[i][a]->ProjectionY(Form("hN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
485 }
486
487 if(fHistPNS[i][a]){
488 (fHistPNS[i][a]->ProjectionY(Form("hPNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
489 (fHistPPS[i][a]->ProjectionY(Form("hPPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
490 (fHistNPS[i][a]->ProjectionY(Form("hNPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
491 (fHistNNS[i][a]->ProjectionY(Form("hNNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
492 (fHistPS[i][a]->ProjectionY(Form("hPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
493 (fHistNS[i][a]->ProjectionY(Form("hNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
494 }
ab34c777 495
496 //printout in text format for delta eta
b7ce3d3f 497 for(Int_t j = 0; j < iCanvas; j++){
ab34c777 498 cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
b7ce3d3f 499
500 if(gbf[i][j][a]){
b7ce3d3f 501 if(a==1){
ab34c777 502 cout<<"Double_t gALICEDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaEta] = {";
503 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
504 cout<<gbf[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
505 cout<<"Double_t gALICEDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaEta] = {";
506 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
507 cout<<gbf[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 508 }
509 else if(a==6){
ab34c777 510 cout<<"Double_t gALICEDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaPhi] = {";
511 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
512 cout<<gbf[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
513 cout<<"Double_t gALICEDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaPhi] = {";
514 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
515 cout<<gbf[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 516 }
517 gbf[i][j][a]->Write();
518 gbf[i][j][a]->Delete();
519 }
520 if(gbfs[i][j][a]){
521 if(a==1){
ab34c777 522 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaEta] = {";
523 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
524 cout<<gbfs[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
525 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaEta] = {";
526 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
527 cout<<gbfs[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 528 }
529 else if(a==6){
ab34c777 530 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaPhi] = {";
531 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
532 cout<<gbfs[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
533 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaPhi] = {";
534 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
535 cout<<gbfs[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 536 }
537 gbfs[i][j][a]->Write();
538 gbfs[i][j][a]->Delete();
539 }
ab34c777 540 cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
541 cout<<endl;
b7ce3d3f 542 }
543 }
544
545 Double_t x,y;
ab34c777 546
547 cout<<"//================================ALICE================================//"<<endl;
b7ce3d3f 548 if(gWM[i]){
ab34c777 549 cout<<"Double_t gWeightedMeanInEtaAlice[nCentralityBins] = {";
550 for(Int_t k = 0; k < gWM[i]->GetN()-1;k++){
b7ce3d3f 551 gWM[i]->GetPoint(k,x,y);
552 cout<<y<<", ";
553 }
ab34c777 554 gWM[i]->GetPoint(k,x,y);
555 cout<<y<<"};"<<endl;
556
557 cout<<"Double_t gWeightedMeanInEtaAliceError[nCentralityBins] = {";
558 for(Int_t k = 0; k < gWM[i]->GetN()-1;k++){
b7ce3d3f 559 cout<<gWM[i]->GetErrorY(k)<<", ";
560 }
ab34c777 561 cout<<gWM[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 562 gWM[i]->Write();
563 }
564 if(gWMS[i]){
ab34c777 565 cout<<"Double_t gShuffledWeightedMeanInEtaAlice[nCentralityBins] = {";
566 for(Int_t k = 0; k < gWMS[i]->GetN()-1;k++){
b7ce3d3f 567 gWMS[i]->GetPoint(k,x,y);
568 cout<<y<<", ";
569 }
ab34c777 570 gWMS[i]->GetPoint(k,x,y);
571 cout<<y<<"};"<<endl;
572
573 cout<<"Double_t gShuffledWeightedMeanInEtaAliceError[nCentralityBins] = {";
574 for(Int_t k = 0; k < gWMS[i]->GetN()-1;k++){
b7ce3d3f 575 cout<<gWMS[i]->GetErrorY(k)<<", ";
576 }
ab34c777 577 cout<<gWMS[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 578 cout<<endl;
579 gWMS[i]->Write();
580 }
581
ab34c777 582 if(gWMP[i]){
583 cout<<"Double_t gWeightedMeanInPhiAlice[nCentralityBins] = {";
584 for(Int_t k = 0; k < gWMP[i]->GetN()-1;k++){
b7ce3d3f 585 gWMP[i]->GetPoint(k,x,y);
586 cout<<y<<", ";
587 }
ab34c777 588 gWMP[i]->GetPoint(k,x,y);
589 cout<<y<<"};"<<endl;
590
591 cout<<"Double_t gWeightedMeanInPhiAliceError[nCentralityBins] = {";
592 for(Int_t k = 0; k < gWMP[i]->GetN()-1;k++){
b7ce3d3f 593 cout<<gWMP[i]->GetErrorY(k)<<", ";
594 }
ab34c777 595 cout<<gWMP[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 596 gWMP[i]->Write();
597 }
598 if(gWMPS[i]){
ab34c777 599 cout<<"Double_t gShuffledWeightedMeanInPhiAlice[nCentralityBins] = {";
600 for(Int_t k = 0; k < gWMPS[i]->GetN()-1;k++){
b7ce3d3f 601 gWMPS[i]->GetPoint(k,x,y);
602 cout<<y<<", ";
603 }
ab34c777 604 gWMPS[i]->GetPoint(k,x,y);
605 cout<<y<<"};"<<endl;
606
607 cout<<"Double_t gShuffledWeightedMeanInPhiAliceError[nCentralityBins] = {";
608 for(Int_t k = 0; k < gWMPS[i]->GetN()-1;k++){
b7ce3d3f 609 cout<<gWMPS[i]->GetErrorY(k)<<", ";
610 }
ab34c777 611 cout<<gWMPS[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 612 cout<<endl;
613 gWMPS[i]->Write();
614 }
615
616 if(ginteg[i]) ginteg[i]->Write();
617 if(gintegS[i]) gintegS[i]->Write();
618 if(gintegP[i]) gintegP[i]->Write();
619 if(gintegPS[i]) gintegPS[i]->Write();
620 }
621 fOut->Close();
622 f->Close();
623 gROOT->Reset();
624}
625
626//____________________________________________________________________//
627void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Double_t &WM, Double_t &WME) {
628
629 //Prints the calculated width of the BF and its error
630 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
631 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
632 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
633 Double_t deltaBalP2 = 0.0, integral = 0.0;
634 Double_t deltaErrorNew = 0.0;
635
636 //Retrieve this variables from Histogram
637 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
638 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
639
640 cout<<"=================================================="<<endl;
641 cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
642 cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
643 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
644 cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
645 }
646 cout<<"=================================================="<<endl;
647 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
648 gSumXi += gHistBalance->GetBinCenter(i);
649 gSumBi += gHistBalance->GetBinContent(i);
650 gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
651 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
652 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
653 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
654 gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
655
656 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
657 integral += fP2Step*gHistBalance->GetBinContent(i);
658 }
659 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
660 deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
661
662 Double_t integralError = TMath::Sqrt(deltaBalP2);
663
664 Double_t delta = gSumBiXi / gSumBi;
665 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
666
667 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
668 cout<<"New error: "<<deltaErrorNew<<endl;
669 cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
670 cout<<"=================================================="<<endl;
671
672 WM = delta;
673 WME = deltaError;
674}
675
676//____________________________________________________________________//
677void GetIntegral(TH1D *gHistBalance, Double_t &integ, Double_t &integE) {
678
679 //always start with 1
680 Int_t fStartBin = 1;
681
682 //Prints the calculated width of the BF and its error
683 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
684 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
685 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
686 Double_t deltaBalP2 = 0.0, integral = 0.0;
687 Double_t deltaErrorNew = 0.0;
688
689 //Retrieve this variables from Histogram
690 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
691 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
692
693 cout<<"=================================================="<<endl;
694 cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
695 cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
696 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
697 cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
698 }
699 cout<<"=================================================="<<endl;
700 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
701 gSumXi += gHistBalance->GetBinCenter(i);
702 gSumBi += gHistBalance->GetBinContent(i);
703 gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
704 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
705 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
706 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
707 gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
708
709 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
710 integral += fP2Step*gHistBalance->GetBinContent(i);
711 }
712 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
713 deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
714
715 Double_t integralError = TMath::Sqrt(deltaBalP2);
716
717 Double_t delta = gSumBiXi / gSumBi;
718 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
719
720 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
721 cout<<"New error: "<<deltaErrorNew<<endl;
722 cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
723 cout<<"=================================================="<<endl;
724
725 integ = integral;
726 integE = integralError;
727}
728
729//___________________________________________________________//
730// NOT USED any more
731//___________________________________________________________//
732void mergeOutput(const char* outputDir) {
733 //Function to merge the output of the sub-jobs
734 //Create a BF object
735 AliBalance *bf = new AliBalance();
736
737 //connect to AliEn's API services
738 TGrid::Connect("alien://");
739
740 //Getting the output dir from the env. variable
741 //(JDL field: JDLVariables={"OutputDir"};)
742 TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");
743
744 Int_t nEntries = result->GetEntries();
745
746 TString alienUrl;
747 TDirectoryFile *dirSubJob;
748
749 TString gCutName[4] = {"Total","Offline trigger",
750 "Vertex","Analyzed"};
751 TH1F *fHistEventStats = new TH1F("fHistEventStats",
752 "Event statistics;;N_{events}",
753 4,0.5,4.5);
754 for(Int_t i = 1; i <= 4; i++)
755 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
756
757 AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");
758 for(Int_t i = 0; i < nEntries; i++) {
759 alienUrl = result->GetKey(i,"turl");
760 alienUrl += "#AnalysisResults.root";
761 Printf("Opening file: %s",alienUrl.Data());
762 TFile *file = TFile::Open(alienUrl.Data());
763 dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWGCFEbyE.outputBalanceFunctionAnalysis.root"));
764
765 //merge BF
766 AliBalance *bfSubJob = dynamic_cast<AliBalance *>(dirSubJob->Get("AliBalance"));
767 //bfSubJob->PrintResults();
768 bf->Merge(bfSubJob);
769 //delete bfSubJob;
770
771 //merge event stats
772 TList *listSubJob = dynamic_cast<TList *>(dirSubJob->Get("listQA"));
773 fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));
774
775 bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));
776 delete listSubJob;
777 }
778
779 //Create the output file
780 TString outputFile = "AnalysisResults.Merged.root";
781 TFile *foutput = TFile::Open(outputFile.Data(),"recreate");
782 TDirectoryFile *dirOutput = new TDirectoryFile();
783 dirOutput->SetName("PWGCFEbyE.outputBalanceFunctionAnalysis.root");
784 //dirOutput->cd();
785 dirOutput->Add(bf);
786 TList *list = new TList();
787 list->SetName("listQA");
788 list->Add(fHistEventStats);
789 list->Add(bfTrackCuts);
790 dirOutput->Add(list);
791 dirOutput->Write();
792 bf->Write();
793 list->Write();
794 foutput->Close();
795
796 //cout<<alienUrl.Data()<<endl;
797}