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