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