]> 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
.so cleanup: removed from gSystem->Load()
[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
47e040b5 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, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE, Bool_t correctWithMixed = kFALSE) {
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
4070f709 16 gSystem->Load("libProofPlayer");
17 gSystem->Load("libANALYSIS");
18 gSystem->Load("libANALYSISalice");
19 gSystem->Load("libEventMixing");
20 gSystem->Load("libPWGCFebye");
b7ce3d3f 21
22 //Draw BF
47e040b5 23 drawBF(bHistos,inFile, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,centEst, "", etaWindow,etaBins, correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed);
b7ce3d3f 24
25
26 //Merge the output
27 //mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");
28}
29
30//___________________________________________________________//
47e040b5 31void 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, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE, Bool_t correctWithMixed = kFALSE) {
b7ce3d3f 32 //Function to draw the BF objects and write them into the output file
33
34 Int_t maximumCanvases = 13;
35 Int_t iCanvas = 0;
36 Int_t iList = -1;
6de53600 37 TCanvas *cQACuts[13][10];
b7ce3d3f 38 TCanvas *cQA[13][10];
39 TCanvas *cQAV0M = new TCanvas("cQAV0M","V0M multiplicities");
40 cQAV0M->Divide(4,3);
41 TCanvas *cQARef = new TCanvas("cQARef","reference track multiplicities");
42 cQARef->Divide(4,3);
43 TCanvas *cBF[13][10];
44 TCanvas *cBFS[13][10];
45
47e040b5 46 // only one correction method allowed (correctWithEfficiency/Acceptance OR divide by event mixing)
47 if((correctWithEfficiency||correctWithAcceptanceOnly) && correctWithMixed){
48 Printf("only one correction method allowed (correctWithEfficiency/Acceptance OR divide by event mixing)");
49 return;
50 }
51
52 // if divide by event mixing, first get all files for all centralities
53 // here the event mixing file is more or less hard coded (has to be changed if somebody else wants to use this)
54 TH1D *hMixedEta[nrOfCentralities][4];
55 TH1D *hMixedPhi[nrOfCentralities][4];
56
57 if(correctWithMixed){
58 TFile* fMixed[nrOfCentralities];
59
60 for (Int_t iCent = 0; iCent < nrOfCentralities ; iCent++){
61 fMixed[iCent] = TFile::Open(Form("/Users/physics/ALICE/balance/BF_EM_new/gridOut/Centrality%.0f_%.0f_oldMethod/Histograms_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_EventMixing_%.0f-%.0f_AnalysisResults.root",centralityArray[iCent],centralityArray[iCent+1],fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,centralityArray[iCent],centralityArray[iCent+1]),"READ");
62
63 if(!fMixed[iCent]){
64 Printf("Event Mixing file not available!");
65 return;
66 }
67
68 // at the time of production the event mixing had no centrality info (all in most central bin, centarlity is selected before)
69 hMixedEta[iCent][0] = (TH1D*)fMixed[iCent]->Get(Form("hPN_eta_0"));
70 hMixedPhi[iCent][0] = (TH1D*)fMixed[iCent]->Get(Form("hPN_phi_0"));
9ff3a85f 71 hMixedEta[iCent][1] = (TH1D*)((TH1D*)fMixed[iCent]->Get(Form("hPN_eta_0")))->Clone("hNP_eta_0");
72 hMixedPhi[iCent][1] = (TH1D*)((TH1D*)fMixed[iCent]->Get(Form("hPN_phi_0")))->Clone("hNP_phi_0");
47e040b5 73 hMixedEta[iCent][2] = (TH1D*)fMixed[iCent]->Get(Form("hNN_eta_0"));
74 hMixedPhi[iCent][2] = (TH1D*)fMixed[iCent]->Get(Form("hNN_phi_0"));
75 hMixedEta[iCent][3] = (TH1D*)fMixed[iCent]->Get(Form("hPP_eta_0"));
76 hMixedPhi[iCent][3] = (TH1D*)fMixed[iCent]->Get(Form("hPP_phi_0"));
9ff3a85f 77
78 // because NP = PN here
79 hMixedEta[iCent][0]->Scale(0.5);
80 hMixedPhi[iCent][0]->Scale(0.5);
81 hMixedEta[iCent][1]->Scale(0.5);
82 hMixedPhi[iCent][1]->Scale(0.5);
83
84 // normalize the distributions to 1:
85 // - in the integral for dphi (for averaging over sector structure)
86 // - in the maximum for deta
87 hMixedPhi[iCent][0]->Scale(1./(Double_t)hMixedPhi[iCent][0]->Integral(hMixedPhi[iCent][0]->FindBin(0),hMixedPhi[iCent][0]->FindBin(180))*(Double_t)(hMixedPhi[iCent][0]->FindBin(180)-hMixedPhi[iCent][0]->FindBin(0)+1));
88 hMixedPhi[iCent][1]->Scale(1./(Double_t)hMixedPhi[iCent][1]->Integral(hMixedPhi[iCent][1]->FindBin(0),hMixedPhi[iCent][1]->FindBin(180))*(Double_t)(hMixedPhi[iCent][1]->FindBin(180)-hMixedPhi[iCent][1]->FindBin(0)+1));
89 hMixedPhi[iCent][2]->Scale(1./(Double_t)hMixedPhi[iCent][2]->Integral(hMixedPhi[iCent][2]->FindBin(0),hMixedPhi[iCent][2]->FindBin(180))*(Double_t)(hMixedPhi[iCent][2]->FindBin(180)-hMixedPhi[iCent][2]->FindBin(0)+1));
90 hMixedPhi[iCent][3]->Scale(1./(Double_t)hMixedPhi[iCent][3]->Integral(hMixedPhi[iCent][3]->FindBin(0),hMixedPhi[iCent][3]->FindBin(180))*(Double_t)(hMixedPhi[iCent][3]->FindBin(180)-hMixedPhi[iCent][3]->FindBin(0)+1));
91
92 hMixedEta[iCent][0]->Scale(0.8 * 20. /(Double_t)hMixedEta[iCent][0]->Integral(hMixedEta[iCent][0]->FindBin(0),hMixedEta[iCent][0]->FindBin(1.6)));
93 hMixedEta[iCent][1]->Scale(0.8 * 20. /(Double_t)hMixedEta[iCent][1]->Integral(hMixedEta[iCent][1]->FindBin(0),hMixedEta[iCent][1]->FindBin(1.6)));
94 hMixedEta[iCent][2]->Scale(0.8 * 20. /(Double_t)hMixedEta[iCent][2]->Integral(hMixedEta[iCent][2]->FindBin(0),hMixedEta[iCent][2]->FindBin(1.6)));
95 hMixedEta[iCent][3]->Scale(0.8 * 20. /(Double_t)hMixedEta[iCent][3]->Integral(hMixedEta[iCent][3]->FindBin(0),hMixedEta[iCent][3]->FindBin(1.6)));
96
47e040b5 97
98 }
99 }
100
b7ce3d3f 101 // get the file
102 TFile *f = TFile::Open(inFile.Data());
103 if(!f) {
104 Printf("File not found!!!");
105 break;
106 }
107
108 // get the BF output directory
109 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis"));
110 if(!dir) {
111 Printf("Output directory not found!!!");
112 break;
113 }
114
115 // loop over all lists and plot the BF and QA
116 TList *list = NULL;
117 TString listName;
118 TIter nextkey( dir->GetListOfKeys() );
119 TKey *key;
120
121 AliBalance *bf[13][7];
122 AliBalance *bfs[13][7];
123 TH1D *gbf[13][10][7];
124 TH1D *gbfs[13][10][7];
125
126 for(Int_t i = 0; i < 13; i++){
127 for(Int_t j = 0; j < 10; j++){
128 for(Int_t k = 0; k < 7; k++){
129 gbf[i][j][k] = NULL;
130 gbfs[i][j][k] = NULL;
131 }
132 }
133 }
134
135 TH2D *fHistP[13][7]; //N+
136 TH2D *fHistN[13][7]; //N-
137 TH2D *fHistPN[13][7]; //N+-
138 TH2D *fHistNP[13][7]; //N-+
139 TH2D *fHistPP[13][7]; //N++
140 TH2D *fHistNN[13][7]; //N--
141
142 TH2D *fHistPS[13][7]; //N+
143 TH2D *fHistNS[13][7]; //N-
144 TH2D *fHistPNS[13][7]; //N+-
145 TH2D *fHistNPS[13][7]; //N-+
146 TH2D *fHistPPS[13][7]; //N++
147 TH2D *fHistNNS[13][7]; //N--
148
149 Double_t WM[13][10]; // weighted mean for eta (recalculated from fStartBin)
150 Double_t WME[13][10]; // error
151 Double_t WMS[13][10]; // weighted mean for eta (recalculated from fStartBin) (shuffled)
152 Double_t WMSE[13][10]; // error (shuffled)
153
ab34c777 154 Double_t WMP[13][10]; // weighted mean for phi (recalculated from fStartBinPhi)
b7ce3d3f 155 Double_t WMPE[13][10]; // error
ab34c777 156 Double_t WMPS[13][10]; // weighted mean for phi (recalculated from fStartBin) (shuffled)
b7ce3d3f 157 Double_t WMPSE[13][10]; // error (shuffled)
158
159 Double_t integ[13][10]; // integral for eta (calculated from bin 1)
160 Double_t integE[13][10]; // error
161 Double_t integS[13][10]; // integral for eta (calculated from bin 1) (shuffled)
162 Double_t integSE[13][10]; // error (shuffled)
163
164 Double_t integP[13][10]; // integral for phi (calculated from bin 1)
165 Double_t integPE[13][10]; // error
166 Double_t integPS[13][10]; // integral for phi (calculated from bin 1) (shuffled)
167 Double_t integPSE[13][10]; // error (shuffled)
168
169
170 while ( (key = (TKey*)nextkey())) {
171
172 list = (TList*)key->ReadObj();
173 listName = TString(list->GetName());
174
47e040b5 175 //cout<<"Processing list "<<listName<<endl;
b7ce3d3f 176
177
178 // ----------------------------------------------------
179 // plot QA histograms
180 if(listName.Contains("QA")){
181
182 if(iList<13) iList++;
183 else{
184 cerr<<"TOO MANY LISTS!!!"<<endl;
185 return;
186 }
187
188 cQA[iList][iCanvas] = new TCanvas(listName,listName);
189 cQA[iList][iCanvas]->Divide(4,3);
190
191 cQA[iList][iCanvas]->cd(1);
192 TH1F* histEventStats = (TH1F*)list->FindObject("fHistEventStats");
193 if(histEventStats){
194 histEventStats->SetFillColor(9);
195 histEventStats->Draw();
196 }
197
198 cQA[iList][iCanvas]->cd(2);
199 TH1F* histTriggerStats = (TH1F*)list->FindObject("fHistTriggerStats");
200 if(histTriggerStats){
201 histTriggerStats->SetFillColor(9);
202 histTriggerStats->Draw();
203 }
204
205 cQA[iList][iCanvas]->cd(3);
206 TH1F* histTrackStats = (TH1F*)list->FindObject("fHistTrackStats");
207 if(histTrackStats){
208 histTrackStats->SetFillColor(9);
209 histTrackStats->Draw();
210 }
211
212 cQA[iList][iCanvas]->cd(4);
213 TH1F* histCentStats = (TH1F*)list->FindObject("fHistCentStats");
214 if(histCentStats){
215 histCentStats->SetFillColor(9);
216 histCentStats->Draw("colz");
217 }
218
219
220
221 cQA[iList][iCanvas]->cd(5);
222 TH1F* histVx = (TH1F*)list->FindObject("fHistVx");
223 if(histVx){
224 histVx->SetFillColor(9);
225 histVx->Draw();
226 }
227 cQA[iList][iCanvas]->cd(6);
228 TH1F* histVy = (TH1F*)list->FindObject("fHistVy");
229 if(histVy){
230 histVy->SetFillColor(9);
231 histVy->Draw();
232 }
233
234 cQA[iList][iCanvas]->cd(7);
235 TH1F* histVz = (TH1F*)list->FindObject("fHistVz");
236 if(histVz){
237 histVz->SetFillColor(9);
238 histVz->Draw();
239 }
240
241 cQA[iList][iCanvas]->cd(8);
242 cQA[iList][iCanvas]->cd(8)->SetLogz();
243 TH2F* histDCA = (TH2F*)list->FindObject("fHistDCA");
244 if(histDCA) histDCA->Draw("colz");
245
246
247 cQA[iList][iCanvas]->cd(9);
248 cQA[iList][iCanvas]->cd(9)->SetLogz();
249 TH2F* histClus = (TH2F*)list->FindObject("fHistClus");
250 if(histClus) histClus->Draw("colz");
251
252 cQA[iList][iCanvas]->cd(10);
253 TH1F* histPt = (TH1F*)list->FindObject("fHistPt");
254 if(histPt){
255 histPt->SetFillColor(9);
256 histPt->Draw();
257 }
258
259 cQA[iList][iCanvas]->cd(11);
260 TH1F* histEta = (TH1F*)list->FindObject("fHistEta");
261 if(histEta){
262 histEta->SetFillColor(9);
263 histEta->Draw();
264 }
265
266 cQA[iList][iCanvas]->cd(12);
267 TH1F* histPhi = (TH1F*)list->FindObject("fHistPhi");
268 if(histPhi){
269 histPhi->SetFillColor(9);
270 histPhi->Draw();
271 }
272
273 // centrality estimator QA
274 cQAV0M->cd(iCanvas+1);
275 cQAV0M->cd(iCanvas+1)->SetLogz();
276 TH1F* histV0M = (TH1F*)list->FindObject("fHistV0M");
277 if(histV0M){
278 histV0M->Draw("colz");
279 }
280
281 cQARef->cd(iCanvas+1);
282 cQARef->cd(iCanvas+1)->SetLogz();
283 TH1F* histRef = (TH1F*)list->FindObject("fHistRefTracks");
284 if(histRef){
285 histRef->Draw("colz");
286 }
287 }
6de53600 288
289
290 cQACuts[iList][iCanvas] = new TCanvas(Form("%sCuts",listName.Data()),Form("%sCuts",listName.Data()));
291 cQACuts[iList][iCanvas]->Divide(2,2);
292
293 cQACuts[iList][iCanvas]->cd(1);
294 TH2D* histHBTcuts = (TH2D*)list->FindObject("fHistHBTafter");
295 if(histHBTcuts){
296 histHBTcuts->DrawCopy("colz");
297 }
298
299 cQACuts[iList][iCanvas]->cd(2);
300 TH2D* histHBTall = (TH2D*)list->FindObject("fHistHBTbefore");
301 if(histHBTall){
302 histHBTcuts->Divide(histHBTall);
303 histHBTcuts->SetMaximum(1.5);
304 histHBTcuts->SetMinimum(0.5);
305 histHBTcuts->DrawCopy("colz");
306 }
307
308 cQACuts[iList][iCanvas]->cd(3);
309 TH2D* histConversioncuts = (TH2D*)list->FindObject("fHistConversionafter");
310 if(histConversioncuts){
311 histConversioncuts->DrawCopy("colz");
312 }
313
314 cQACuts[iList][iCanvas]->cd(4);
315 TH2D* histConversionall = (TH2D*)list->FindObject("fHistConversionbefore");
316 if(histConversionall){
317 histConversioncuts->Divide(histConversionall);
318 histConversioncuts->SetMaximum(1.5);
319 histConversioncuts->SetMinimum(0.5);
320 histConversioncuts->DrawCopy("colz");
321 }
322
b7ce3d3f 323 // ----------------------------------------------------
324
325 // ----------------------------------------------------
326 // calculate and plot BF
ab34c777 327 if(listName.Contains("BF_")&&listName.Contains(centEst.Data())){
b7ce3d3f 328
329 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
330
331 cBF[iList][iCanvas] = new TCanvas(Form("cBF_%d_%d",iList,iCanvas),Form("cBF_%d_%d",iList,iCanvas));
332 cBF[iList][iCanvas]->Divide(3,3);
333
334 }
335
336
337 for(Int_t a = 0; a < 7; a++){
338
47e040b5 339 //cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;
b7ce3d3f 340
341 // create the BF object
342 bf[iList][a] = new AliBalance();
343
344 fHistP[iList][a] = (TH2D*)list->FindObject(Form("fHistP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
345 fHistN[iList][a] = (TH2D*)list->FindObject(Form("fHistN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
346 fHistPP[iList][a] = (TH2D*)list->FindObject(Form("fHistPP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
347 fHistPN[iList][a] = (TH2D*)list->FindObject(Form("fHistPN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
348 fHistNP[iList][a] = (TH2D*)list->FindObject(Form("fHistNP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
349 fHistNN[iList][a] = (TH2D*)list->FindObject(Form("fHistNN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
350
351 // rebin histograms (be careful with divider!)
352 if(a==6){
ab34c777 353 fHistP[iList][a]->RebinY(fRebinPhi);
354 fHistN[iList][a]->RebinY(fRebinPhi);
355 fHistPP[iList][a]->RebinY(fRebinPhi);
356 fHistPN[iList][a]->RebinY(fRebinPhi);
357 fHistNP[iList][a]->RebinY(fRebinPhi);
358 fHistNN[iList][a]->RebinY(fRebinPhi);
b7ce3d3f 359 }
360 else{
361 fHistP[iList][a]->RebinY(fRebin);
362 fHistN[iList][a]->RebinY(fRebin);
363 fHistPP[iList][a]->RebinY(fRebin);
364 fHistPN[iList][a]->RebinY(fRebin);
365 fHistNP[iList][a]->RebinY(fRebin);
366 fHistNN[iList][a]->RebinY(fRebin);
367 }
368
369 fHistP[iList][a]->SetName(Form("%s_%d_%d",fHistP[iList][a]->GetName(),iList,iCanvas));
370 fHistN[iList][a]->SetName(Form("%s_%d_%d",fHistN[iList][a]->GetName(),iList,iCanvas));
371 fHistPP[iList][a]->SetName(Form("%s_%d_%d",fHistPP[iList][a]->GetName(),iList,iCanvas));
372 fHistPN[iList][a]->SetName(Form("%s_%d_%d",fHistPN[iList][a]->GetName(),iList,iCanvas));
373 fHistNP[iList][a]->SetName(Form("%s_%d_%d",fHistNP[iList][a]->GetName(),iList,iCanvas));
374 fHistNN[iList][a]->SetName(Form("%s_%d_%d",fHistNN[iList][a]->GetName(),iList,iCanvas));
375
376 // set histograms in AliBalance object
377 bf[iList][a]->SetHistNp(a, fHistP[iList][a]);
378 bf[iList][a]->SetHistNn(a, fHistN[iList][a]);
379 bf[iList][a]->SetHistNpp(a, fHistPP[iList][a]);
380 bf[iList][a]->SetHistNpn(a, fHistPN[iList][a]);
381 bf[iList][a]->SetHistNnp(a, fHistNP[iList][a]);
382 bf[iList][a]->SetHistNnn(a, fHistNN[iList][a]);
383
384 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
47e040b5 385
386 if(correctWithMixed && a == 1) gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow,correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedEta[iCanvas]);
387 else if(correctWithMixed && a == 6) gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow,correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedPhi[iCanvas]);
388 else gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow,correctWithEfficiency,correctWithAcceptanceOnly);
b7ce3d3f 389
b7ce3d3f 390 gbf[iList][iCanvas][a]->SetName(Form("BF_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
391
392 cBF[iList][iCanvas]->cd(a+1);
393 gbf[iList][iCanvas][a]->SetMarkerStyle(20);
394
395 if(!bHistos){
396 gbf[iList][iCanvas][a]->DrawCopy("AP");
397 if(a==1){
ce2aff27 398 GetWeightedMean(gbf[iList][iCanvas][a],fStartBinBFWidth,etaBins,WM[iList][iCanvas],WME[iList][iCanvas]); // for eta recalculate width
b7ce3d3f 399 GetIntegral(gbf[iList][iCanvas][a],integ[iList][iCanvas],integE[iList][iCanvas]);
400 }
401 else if(a==6){
ce2aff27 402 GetWeightedMean(gbf[iList][iCanvas][a],fStartBinBFWidthPhi,-1,WMP[iList][iCanvas],WMPE[iList][iCanvas]); // for phi calculate width
b7ce3d3f 403 GetIntegral(gbf[iList][iCanvas][a],integP[iList][iCanvas],integPE[iList][iCanvas]);
404 }
405 }
406 else{
407 fHistPN[iList][a]->SetLineColor(2);
6de53600 408 TH1D *fHistTmp = (TH1D*)fHistPN[iList][a]->ProjectionY(Form("pn%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1]);
409 fHistTmp->Scale(0.5);
410 fHistTmp->DrawCopy();
b7ce3d3f 411 fHistPP[iList][a]->SetLineColor(1);
6de53600 412 fHistPP[iList][a]->ProjectionY(Form("pp%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
b7ce3d3f 413 fHistNP[iList][a]->SetLineColor(4);
6de53600 414 fHistNP[iList][a]->ProjectionY(Form("np%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
b7ce3d3f 415 fHistNN[iList][a]->SetLineColor(8);
6de53600 416 fHistNN[iList][a]->ProjectionY(Form("nn%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
b7ce3d3f 417 }
418 }
419 }
6de53600 420
421 if(bHistos){
422 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
423 cBF[iList][iCanvas]->cd(8);
424 fHistP[iList][1]->ProjectionY(Form("p%d",1),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy();
425 fHistN[iList][1]->ProjectionY(Form("n%d",1),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
426 cBF[iList][iCanvas]->cd(9);
427 fHistP[iList][6]->ProjectionY(Form("p%d",6),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy();
428 fHistN[iList][6]->ProjectionY(Form("n%d",6),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
429 }
430 }
b7ce3d3f 431 }
432 // ----------------------------------------------------
433
434 // ----------------------------------------------------
435 // calculate and plot BF (shuffled)
ab34c777 436 if(listName.Contains("BFShuffled")&&listName.Contains(centEst.Data())&&listName.Contains(extraString.Data())){
b7ce3d3f 437
438 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
439
440 cBFS[iList][iCanvas] = new TCanvas(Form("Shuffled_%d_%d",iList,iCanvas),Form("Shuffled_%d_%d",iList,iCanvas));
441 cBFS[iList][iCanvas]->Divide(3,3);
442
443 }
444
445 for(Int_t a = 0; a < 7; a++){
446
447 // create the BF object
448 bfs[iList][a] = new AliBalance();
449
450 fHistPS[iList][a] = (TH2D*)list->FindObject(Form("fHistP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
451 fHistNS[iList][a] = (TH2D*)list->FindObject(Form("fHistN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
452 fHistPPS[iList][a] = (TH2D*)list->FindObject(Form("fHistPP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
453 fHistPNS[iList][a] = (TH2D*)list->FindObject(Form("fHistPN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
454 fHistNPS[iList][a] = (TH2D*)list->FindObject(Form("fHistNP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
455 fHistNNS[iList][a] = (TH2D*)list->FindObject(Form("fHistNN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
456
457 // rebin histograms (be careful with divider!)
458 if(a==6){
ab34c777 459 fHistPS[iList][a]->RebinY(fRebinPhi);
460 fHistNS[iList][a]->RebinY(fRebinPhi);
461 fHistPPS[iList][a]->RebinY(fRebinPhi);
462 fHistPNS[iList][a]->RebinY(fRebinPhi);
463 fHistNPS[iList][a]->RebinY(fRebinPhi);
464 fHistNNS[iList][a]->RebinY(fRebinPhi);
b7ce3d3f 465 }
466 else{
467 fHistPS[iList][a]->RebinY(fRebin);
468 fHistNS[iList][a]->RebinY(fRebin);
469 fHistPPS[iList][a]->RebinY(fRebin);
470 fHistPNS[iList][a]->RebinY(fRebin);
471 fHistNPS[iList][a]->RebinY(fRebin);
472 fHistNNS[iList][a]->RebinY(fRebin);
473 }
474
475 fHistPS[iList][a]->SetName(Form("%s_%d_%d",fHistPS[iList][a]->GetName(),iList,iCanvas));
476 fHistNS[iList][a]->SetName(Form("%s_%d_%d",fHistNS[iList][a]->GetName(),iList,iCanvas));
477 fHistPPS[iList][a]->SetName(Form("%s_%d_%d",fHistPPS[iList][a]->GetName(),iList,iCanvas));
478 fHistPNS[iList][a]->SetName(Form("%s_%d_%d",fHistPNS[iList][a]->GetName(),iList,iCanvas));
479 fHistNPS[iList][a]->SetName(Form("%s_%d_%d",fHistNPS[iList][a]->GetName(),iList,iCanvas));
480 fHistNNS[iList][a]->SetName(Form("%s_%d_%d",fHistNNS[iList][a]->GetName(),iList,iCanvas));
481
482 // set histograms in AliBalance object
483 bfs[iList][a]->SetHistNp(a, fHistPS[iList][a]);
484 bfs[iList][a]->SetHistNn(a, fHistNS[iList][a]);
485 bfs[iList][a]->SetHistNpp(a, fHistPPS[iList][a]);
486 bfs[iList][a]->SetHistNpn(a, fHistPNS[iList][a]);
487 bfs[iList][a]->SetHistNnp(a, fHistNPS[iList][a]);
488 bfs[iList][a]->SetHistNnn(a, fHistNNS[iList][a]);
489
490 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
491
47e040b5 492 if(correctWithMixed && a == 1) gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow, correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedEta[iCanvas]);
493 else if(correctWithMixed && a == 6) gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow, correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedPhi[iCanvas]);
494 else gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow, correctWithEfficiency,correctWithAcceptanceOnly);
b7ce3d3f 495 gbfs[iList][iCanvas][a]->SetName(Form("BFS_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
496
497 cBFS[iList][iCanvas]->cd(a+1);
498 gbfs[iList][iCanvas][a]->SetMarkerStyle(20);
499 if(!bHistos){
500 gbfs[iList][iCanvas][a]->DrawCopy("AP");
501 if(a==1){
ce2aff27 502 GetWeightedMean(gbfs[iList][iCanvas][a],fStartBinBFWidth,etaBins,WMS[iList][iCanvas],WMSE[iList][iCanvas]);
b7ce3d3f 503 GetIntegral(gbfs[iList][iCanvas][a],integS[iList][iCanvas],integSE[iList][iCanvas]);
504 }
505 else if(a==6){
ce2aff27 506 GetWeightedMean(gbfs[iList][iCanvas][a],fStartBinBFWidthPhi,-1,WMPS[iList][iCanvas],WMPSE[iList][iCanvas]); // for phi calculate width
b7ce3d3f 507 GetIntegral(gbfs[iList][iCanvas][a],integPS[iList][iCanvas],integPSE[iList][iCanvas]);
508 }
509 }
510 else{
511 fHistPNS[iList][a]->SetLineColor(2);
6de53600 512 TH1D *fHistTmp = (TH1D*)fHistPN[iList][a]->ProjectionY(Form("pns%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1]);
513 fHistTmp->Scale(0.5);
514 fHistTmp->DrawCopy();
b7ce3d3f 515 fHistPPS[iList][a]->SetLineColor(1);
6de53600 516 fHistPPS[iList][a]->ProjectionY(Form("pps%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
b7ce3d3f 517 fHistNPS[iList][a]->SetLineColor(4);
6de53600 518 fHistNPS[iList][a]->ProjectionY(Form("nps%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
b7ce3d3f 519 fHistNNS[iList][a]->SetLineColor(8);
6de53600 520 fHistNNS[iList][a]->ProjectionY(Form("nns%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
b7ce3d3f 521 }
522 }
523 }
6de53600 524
525 if(bHistos){
526 for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
527 cBFS[iList][iCanvas]->cd(8);
528 fHistPS[iList][1]->ProjectionY(Form("pS%d",1),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy();
529 fHistNS[iList][1]->ProjectionY(Form("nS%d",1),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
530 cBFS[iList][iCanvas]->cd(9);
531 fHistPS[iList][6]->ProjectionY(Form("pS%d",6),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy();
532 fHistNS[iList][6]->ProjectionY(Form("nS%d",6),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
533 }
534 }
b7ce3d3f 535 }
536 // ----------------------------------------------------
537 }
538
539 // for BF calculation create also graphs with weighted mean for eta
540 if(iCanvas == 0) return;
541
542 TGraphErrors *gWM[13];
543 TGraphErrors *gWMS[13];
544 TGraphErrors *gWMP[13];
545 TGraphErrors *gWMPS[13];
546 TGraphErrors *ginteg[13];
547 TGraphErrors *gintegS[13];
548 TGraphErrors *gintegP[13];
549 TGraphErrors *gintegPS[13];
550 for(Int_t i = 0; i < 13; i++){
551 gWM[i] = NULL;
552 gWMS[i] = NULL;
553 ginteg[i] = NULL;
554 gintegS[i] = NULL;
555 gWMP[i] = NULL;
556 gWMPS[i] = NULL;
557 gintegP[i] = NULL;
558 gintegPS[i] = NULL;
559 }
560
561
562 if(!bHistos){
563 for(Int_t i = 0; i < iList+1; i++){
564 gWM[i] = new TGraphErrors(iCanvas,cent,WM[i],centE,WME[i]);
036e5409 565 if(iList==0) gWM[i]->SetName("gCentrality");
b7ce3d3f 566 else gWM[i]->SetName(Form("gCentrality_%d",i));
567 gWMS[i] = new TGraphErrors(iCanvas,cent,WMS[i],centE,WMSE[i]);
036e5409 568 if(iList==0) gWMS[i]->SetName("gCentralityS");
b7ce3d3f 569 else gWMS[i]->SetName(Form("gCentralityS_%d",i));
570
571 gWMP[i] = new TGraphErrors(iCanvas,cent,WMP[i],centE,WMPE[i]);
036e5409 572 if(iList==0) gWMP[i]->SetName("gCentralityPhi");
b7ce3d3f 573 else gWMP[i]->SetName(Form("gCentralityPhi_%d",i));
574 gWMPS[i] = new TGraphErrors(iCanvas,cent,WMPS[i],centE,WMPSE[i]);
036e5409 575 if(iList==0) gWMPS[i]->SetName("gCentralityPhiS");
b7ce3d3f 576 else gWMPS[i]->SetName(Form("gCentralityPhiS_%d",i));
577
578 ginteg[i] = new TGraphErrors(iCanvas,cent,integ[i],centE,integE[i]);
036e5409 579 if(iList==0) ginteg[i]->SetName("gIntegral");
b7ce3d3f 580 else ginteg[i]->SetName(Form("gIntegral_%d",i));
581 gintegS[i] = new TGraphErrors(iCanvas,cent,integS[i],centE,integSE[i]);
036e5409 582 if(iList==0) gintegS[i]->SetName("gIntegralS");
b7ce3d3f 583 else gintegS[i]->SetName(Form("gIntegralS_%d",i));
584
585 gintegP[i] = new TGraphErrors(iCanvas,cent,integP[i],centE,integPE[i]);
036e5409 586 if(iList==0) gintegP[i]->SetName("gIntegraPhil");
b7ce3d3f 587 else gintegP[i]->SetName(Form("gIntegralPhi_%d",i));
588 gintegPS[i] = new TGraphErrors(iCanvas,cent,integPS[i],centE,integPSE[i]);
036e5409 589 if(iList==0) gintegPS[i]->SetName("gIntegralPhiS");
b7ce3d3f 590 else gintegPS[i]->SetName(Form("gIntegralPhiS_%d",i));
591 }
592 }
593
ce2aff27 594 TFile *fOut = NULL;
595 if(etaWindow > 0 || etaBins > -1){
1c8a00af 596 if(correctWithEfficiency){
b5b9a4f9 597 if(correctWithAcceptanceOnly){
598 fOut = TFile::Open(Form("Histograms_AccCorrWithAcceptanceOnly_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
599 }
600 else{
601 fOut = TFile::Open(Form("Histograms_AccCorrWithEfficiency_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
602 }
1c8a00af 603 }
47e040b5 604 else if(correctWithMixed){
605 fOut = TFile::Open(Form("Histograms_CorrEM_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
606 }
1c8a00af 607 else{
608 fOut = TFile::Open(Form("Histograms_AccCorr_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
609 }
ce2aff27 610 }
611 else{
612 fOut = TFile::Open(Form("Histograms_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
613 }
b7ce3d3f 614 fOut->cd();
615 for(Int_t i = 0; i < iList+1; i++){
47e040b5 616 if(gbf[i][0][1]){
ab34c777 617 cout<<"PROCESS LIST "<<i<<" NOW!"<<endl;
b7ce3d3f 618 for(Int_t a = 0; a < 7; a++){
47e040b5 619 if(a==1||a==6){
620 cout<<"PROCESS VARIABLE "<<a<<" NOW!"<<endl;
b7ce3d3f 621
622 if(fHistPN[i][a]){
623 (fHistPN[i][a]->ProjectionY(Form("hPN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
624 (fHistPP[i][a]->ProjectionY(Form("hPP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
625 (fHistNP[i][a]->ProjectionY(Form("hNP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
626 (fHistNN[i][a]->ProjectionY(Form("hNN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
627 (fHistP[i][a]->ProjectionY(Form("hP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
628 (fHistN[i][a]->ProjectionY(Form("hN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
629 }
630
631 if(fHistPNS[i][a]){
632 (fHistPNS[i][a]->ProjectionY(Form("hPNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
633 (fHistPPS[i][a]->ProjectionY(Form("hPPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
634 (fHistNPS[i][a]->ProjectionY(Form("hNPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
635 (fHistNNS[i][a]->ProjectionY(Form("hNNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
636 (fHistPS[i][a]->ProjectionY(Form("hPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
637 (fHistNS[i][a]->ProjectionY(Form("hNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
638 }
ab34c777 639
640 //printout in text format for delta eta
b7ce3d3f 641 for(Int_t j = 0; j < iCanvas; j++){
ab34c777 642 cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
b7ce3d3f 643
644 if(gbf[i][j][a]){
b7ce3d3f 645 if(a==1){
ab34c777 646 cout<<"Double_t gALICEDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaEta] = {";
647 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
648 cout<<gbf[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
649 cout<<"Double_t gALICEDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaEta] = {";
650 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
651 cout<<gbf[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 652 }
653 else if(a==6){
ab34c777 654 cout<<"Double_t gALICEDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaPhi] = {";
655 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
656 cout<<gbf[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
657 cout<<"Double_t gALICEDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaPhi] = {";
658 for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
659 cout<<gbf[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 660 }
661 gbf[i][j][a]->Write();
662 gbf[i][j][a]->Delete();
663 }
664 if(gbfs[i][j][a]){
665 if(a==1){
ab34c777 666 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaEta] = {";
667 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
668 cout<<gbfs[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
669 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaEta] = {";
670 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
671 cout<<gbfs[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 672 }
673 else if(a==6){
ab34c777 674 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaPhi] = {";
675 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
676 cout<<gbfs[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
677 cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaPhi] = {";
678 for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
679 cout<<gbfs[i][j][a]->GetBinError(k+1)<<"};"<<endl;
b7ce3d3f 680 }
681 gbfs[i][j][a]->Write();
682 gbfs[i][j][a]->Delete();
683 }
ab34c777 684 cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
685 cout<<endl;
b7ce3d3f 686 }
47e040b5 687 }
b7ce3d3f 688 }
689
690 Double_t x,y;
ab34c777 691
692 cout<<"//================================ALICE================================//"<<endl;
b7ce3d3f 693 if(gWM[i]){
ab34c777 694 cout<<"Double_t gWeightedMeanInEtaAlice[nCentralityBins] = {";
695 for(Int_t k = 0; k < gWM[i]->GetN()-1;k++){
b7ce3d3f 696 gWM[i]->GetPoint(k,x,y);
697 cout<<y<<", ";
698 }
ab34c777 699 gWM[i]->GetPoint(k,x,y);
700 cout<<y<<"};"<<endl;
701
702 cout<<"Double_t gWeightedMeanInEtaAliceError[nCentralityBins] = {";
703 for(Int_t k = 0; k < gWM[i]->GetN()-1;k++){
b7ce3d3f 704 cout<<gWM[i]->GetErrorY(k)<<", ";
705 }
ab34c777 706 cout<<gWM[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 707 gWM[i]->Write();
708 }
709 if(gWMS[i]){
ab34c777 710 cout<<"Double_t gShuffledWeightedMeanInEtaAlice[nCentralityBins] = {";
711 for(Int_t k = 0; k < gWMS[i]->GetN()-1;k++){
b7ce3d3f 712 gWMS[i]->GetPoint(k,x,y);
713 cout<<y<<", ";
714 }
ab34c777 715 gWMS[i]->GetPoint(k,x,y);
716 cout<<y<<"};"<<endl;
717
718 cout<<"Double_t gShuffledWeightedMeanInEtaAliceError[nCentralityBins] = {";
719 for(Int_t k = 0; k < gWMS[i]->GetN()-1;k++){
b7ce3d3f 720 cout<<gWMS[i]->GetErrorY(k)<<", ";
721 }
ab34c777 722 cout<<gWMS[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 723 cout<<endl;
724 gWMS[i]->Write();
725 }
726
ab34c777 727 if(gWMP[i]){
728 cout<<"Double_t gWeightedMeanInPhiAlice[nCentralityBins] = {";
729 for(Int_t k = 0; k < gWMP[i]->GetN()-1;k++){
b7ce3d3f 730 gWMP[i]->GetPoint(k,x,y);
731 cout<<y<<", ";
732 }
ab34c777 733 gWMP[i]->GetPoint(k,x,y);
734 cout<<y<<"};"<<endl;
735
736 cout<<"Double_t gWeightedMeanInPhiAliceError[nCentralityBins] = {";
737 for(Int_t k = 0; k < gWMP[i]->GetN()-1;k++){
b7ce3d3f 738 cout<<gWMP[i]->GetErrorY(k)<<", ";
739 }
ab34c777 740 cout<<gWMP[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 741 gWMP[i]->Write();
742 }
743 if(gWMPS[i]){
ab34c777 744 cout<<"Double_t gShuffledWeightedMeanInPhiAlice[nCentralityBins] = {";
745 for(Int_t k = 0; k < gWMPS[i]->GetN()-1;k++){
b7ce3d3f 746 gWMPS[i]->GetPoint(k,x,y);
747 cout<<y<<", ";
748 }
ab34c777 749 gWMPS[i]->GetPoint(k,x,y);
750 cout<<y<<"};"<<endl;
751
752 cout<<"Double_t gShuffledWeightedMeanInPhiAliceError[nCentralityBins] = {";
753 for(Int_t k = 0; k < gWMPS[i]->GetN()-1;k++){
b7ce3d3f 754 cout<<gWMPS[i]->GetErrorY(k)<<", ";
755 }
ab34c777 756 cout<<gWMPS[i]->GetErrorY(k)<<"};"<<endl;
b7ce3d3f 757 cout<<endl;
758 gWMPS[i]->Write();
c418774f 759 }
760 if(ginteg[i]){
761 cout<<"Double_t gIntegralEtaAlice[nCentralityBins] = {";
762 for(Int_t k = 0; k < ginteg[i]->GetN()-1;k++){
763 ginteg[i]->GetPoint(k,x,y);
764 cout<<y<<", ";
765 }
766 ginteg[i]->GetPoint(k,x,y);
767 cout<<y<<"};"<<endl;
768 cout<<endl;
769
770 cout<<"Double_t gIntegralErrorEtaAlice[nCentralityBins] = {";
771 for(Int_t k = 0; k < ginteg[i]->GetN()-1;k++){
772 cout<<ginteg[i]->GetErrorY(k)<<", ";
773 }
774 cout<<ginteg[i]->GetErrorY(k)<<"};"<<endl;
775 cout<<endl;
776
777 ginteg[i]->Write();
778
779 cout<<"Double_t gIntegralPhiAlice[nCentralityBins] = {";
780 for(Int_t k = 0; k < gintegP[i]->GetN()-1;k++){
781 gintegP[i]->GetPoint(k,x,y);
782 cout<<y<<", ";
783 }
784 gintegP[i]->GetPoint(k,x,y);
785 cout<<y<<"};"<<endl;
786 cout<<endl;
787
788 cout<<"Double_t gIntegralErrorPhiAlice[nCentralityBins] = {";
789 for(Int_t k = 0; k < gintegP[i]->GetN()-1;k++){
790 cout<<gintegP[i]->GetErrorY(k)<<", ";
791 }
792 cout<<gintegP[i]->GetErrorY(k)<<"};"<<endl;
793 cout<<endl;
794
795 gintegP[i]->Write();
b7ce3d3f 796
c418774f 797 cout<<"Double_t gIntegralEtaShuffledAlice[nCentralityBins] = {";
798 for(Int_t k = 0; k < gintegS[i]->GetN()-1;k++){
799 gintegS[i]->GetPoint(k,x,y);
800 cout<<y<<", ";
801 }
802 gintegS[i]->GetPoint(k,x,y);
803 cout<<y<<"};"<<endl;
804 cout<<endl;
805
806 cout<<"Double_t gIntegralErrorEtaShuffledAlice[nCentralityBins] = {";
807 for(Int_t k = 0; k < gintegS[i]->GetN()-1;k++){
808 cout<<gintegS[i]->GetErrorY(k)<<", ";
809 }
810 cout<<gintegS[i]->GetErrorY(k)<<"};"<<endl;
811 cout<<endl;
812
813 gintegS[i]->Write();
814
815 cout<<"Double_t gIntegralPhiShuffledAlice[nCentralityBins] = {";
816 for(Int_t k = 0; k < gintegPS[i]->GetN()-1;k++){
3f298f0b 817 gintegPS[i]->GetPoint(k,x,y);
c418774f 818 cout<<y<<", ";
819 }
820 gintegPS[i]->GetPoint(k,x,y);
821 cout<<y<<"};"<<endl;
822 cout<<endl;
823
824 cout<<"Double_t gIntegralErrorPhiShuffledAlice[nCentralityBins] = {";
825 for(Int_t k = 0; k < gintegPS[i]->GetN()-1;k++){
826 cout<<gintegPS[i]->GetErrorY(k)<<", ";
827 }
828 cout<<gintegPS[i]->GetErrorY(k)<<"};"<<endl;
829 cout<<endl;
830
831 gintegPS[i]->Write();
832
833 }
47e040b5 834 }
b7ce3d3f 835 }
836 fOut->Close();
837 f->Close();
838 gROOT->Reset();
839}
840
841//____________________________________________________________________//
ce2aff27 842void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME) {
b7ce3d3f 843
844 //Prints the calculated width of the BF and its error
845 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
846 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
847 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
848 Double_t deltaBalP2 = 0.0, integral = 0.0;
849 Double_t deltaErrorNew = 0.0;
850
851 //Retrieve this variables from Histogram
852 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
ce2aff27 853 if(fStopBin > -1) fNumberOfBins = fStopBin;
b7ce3d3f 854 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
855
47e040b5 856 // cout<<"=================================================="<<endl;
857 // cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
858 // cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
859 // for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
860 // cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
861 // }
862 // cout<<"=================================================="<<endl;
b7ce3d3f 863 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
864 gSumXi += gHistBalance->GetBinCenter(i);
865 gSumBi += gHistBalance->GetBinContent(i);
866 gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
867 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
868 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
869 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
870 gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
871
872 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
873 integral += fP2Step*gHistBalance->GetBinContent(i);
874 }
875 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
876 deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
877
878 Double_t integralError = TMath::Sqrt(deltaBalP2);
879
880 Double_t delta = gSumBiXi / gSumBi;
881 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
882
47e040b5 883 // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
884 // cout<<"New error: "<<deltaErrorNew<<endl;
885 // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
886 // cout<<"=================================================="<<endl;
b7ce3d3f 887
888 WM = delta;
889 WME = deltaError;
890}
891
892//____________________________________________________________________//
893void GetIntegral(TH1D *gHistBalance, Double_t &integ, Double_t &integE) {
894
895 //always start with 1
896 Int_t fStartBin = 1;
897
898 //Prints the calculated width of the BF and its error
899 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
900 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
901 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
902 Double_t deltaBalP2 = 0.0, integral = 0.0;
903 Double_t deltaErrorNew = 0.0;
904
905 //Retrieve this variables from Histogram
906 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
907 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
908
47e040b5 909 // cout<<"=================================================="<<endl;
910 // cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
911 // cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
912 // for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
913 // cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
914 // }
915 // cout<<"=================================================="<<endl;
b7ce3d3f 916 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
917 gSumXi += gHistBalance->GetBinCenter(i);
918 gSumBi += gHistBalance->GetBinContent(i);
919 gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
920 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
921 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
922 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
923 gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
924
925 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
926 integral += fP2Step*gHistBalance->GetBinContent(i);
927 }
928 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
929 deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
930
931 Double_t integralError = TMath::Sqrt(deltaBalP2);
932
933 Double_t delta = gSumBiXi / gSumBi;
934 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
935
47e040b5 936 // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
937 // cout<<"New error: "<<deltaErrorNew<<endl;
938 // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
939 // cout<<"=================================================="<<endl;
b7ce3d3f 940
941 integ = integral;
942 integE = integralError;
943}
944
945//___________________________________________________________//
946// NOT USED any more
947//___________________________________________________________//
948void mergeOutput(const char* outputDir) {
949 //Function to merge the output of the sub-jobs
950 //Create a BF object
951 AliBalance *bf = new AliBalance();
952
953 //connect to AliEn's API services
954 TGrid::Connect("alien://");
955
956 //Getting the output dir from the env. variable
957 //(JDL field: JDLVariables={"OutputDir"};)
958 TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");
959
960 Int_t nEntries = result->GetEntries();
961
962 TString alienUrl;
963 TDirectoryFile *dirSubJob;
964
965 TString gCutName[4] = {"Total","Offline trigger",
966 "Vertex","Analyzed"};
967 TH1F *fHistEventStats = new TH1F("fHistEventStats",
968 "Event statistics;;N_{events}",
969 4,0.5,4.5);
970 for(Int_t i = 1; i <= 4; i++)
971 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
972
973 AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");
974 for(Int_t i = 0; i < nEntries; i++) {
975 alienUrl = result->GetKey(i,"turl");
976 alienUrl += "#AnalysisResults.root";
977 Printf("Opening file: %s",alienUrl.Data());
978 TFile *file = TFile::Open(alienUrl.Data());
979 dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWGCFEbyE.outputBalanceFunctionAnalysis.root"));
980
981 //merge BF
982 AliBalance *bfSubJob = dynamic_cast<AliBalance *>(dirSubJob->Get("AliBalance"));
983 //bfSubJob->PrintResults();
984 bf->Merge(bfSubJob);
985 //delete bfSubJob;
986
987 //merge event stats
988 TList *listSubJob = dynamic_cast<TList *>(dirSubJob->Get("listQA"));
989 fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));
990
991 bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));
992 delete listSubJob;
993 }
994
995 //Create the output file
996 TString outputFile = "AnalysisResults.Merged.root";
997 TFile *foutput = TFile::Open(outputFile.Data(),"recreate");
998 TDirectoryFile *dirOutput = new TDirectoryFile();
999 dirOutput->SetName("PWGCFEbyE.outputBalanceFunctionAnalysis.root");
1000 //dirOutput->cd();
1001 dirOutput->Add(bf);
1002 TList *list = new TList();
1003 list->SetName("listQA");
1004 list->Add(fHistEventStats);
1005 list->Add(bfTrackCuts);
1006 dirOutput->Add(list);
1007 dirOutput->Write();
1008 bf->Write();
1009 list->Write();
1010 foutput->Close();
1011
1012 //cout<<alienUrl.Data()<<endl;
1013}