.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / readBalanceFunction.C
1 const TString gBFAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
2
3 const Int_t nrOfCentralities = 9;
4 const 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)
6 const 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 
7 const 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
9 void 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) {
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");
17   gSystem->Load("libANALYSIS");
18   gSystem->Load("libANALYSISalice");
19   gSystem->Load("libEventMixing");
20   gSystem->Load("libPWGCFebye");
21
22   //Draw BF       
23   drawBF(bHistos,inFile, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,centEst, "",  etaWindow,etaBins, correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed);    
24   
25
26   //Merge the output
27   //mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");
28 }
29
30 //___________________________________________________________//
31 void 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) {
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;
37   TCanvas *cQACuts[13][10];
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   // 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*)((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");
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       // 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     
97       
98     }
99   }
100
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
154   Double_t WMP[13][10];     // weighted mean for phi (recalculated from fStartBinPhi)
155   Double_t WMPE[13][10];    // error
156   Double_t WMPS[13][10];     // weighted mean for phi (recalculated from fStartBin) (shuffled)
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
175     //cout<<"Processing list "<<listName<<endl;
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     }
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
323     // ----------------------------------------------------
324
325     // ----------------------------------------------------
326     // calculate and plot BF 
327     if(listName.Contains("BF_")&&listName.Contains(centEst.Data())){
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
339         //cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;
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){
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);
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++){
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);
389
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){
398               GetWeightedMean(gbf[iList][iCanvas][a],fStartBinBFWidth,etaBins,WM[iList][iCanvas],WME[iList][iCanvas]); // for eta recalculate width 
399               GetIntegral(gbf[iList][iCanvas][a],integ[iList][iCanvas],integE[iList][iCanvas]);
400             }
401             else if(a==6){
402               GetWeightedMean(gbf[iList][iCanvas][a],fStartBinBFWidthPhi,-1,WMP[iList][iCanvas],WMPE[iList][iCanvas]); // for phi calculate width 
403               GetIntegral(gbf[iList][iCanvas][a],integP[iList][iCanvas],integPE[iList][iCanvas]);
404             }
405           }
406           else{
407             fHistPN[iList][a]->SetLineColor(2);
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();
411             fHistPP[iList][a]->SetLineColor(1);
412             fHistPP[iList][a]->ProjectionY(Form("pp%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
413             fHistNP[iList][a]->SetLineColor(4);
414             fHistNP[iList][a]->ProjectionY(Form("np%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
415             fHistNN[iList][a]->SetLineColor(8);
416             fHistNN[iList][a]->ProjectionY(Form("nn%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
417           }
418         }
419       }
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       }
431     }
432     // ----------------------------------------------------
433
434     // ----------------------------------------------------
435     // calculate and plot BF (shuffled)
436     if(listName.Contains("BFShuffled")&&listName.Contains(centEst.Data())&&listName.Contains(extraString.Data())){
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){
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);
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           
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);
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){
502               GetWeightedMean(gbfs[iList][iCanvas][a],fStartBinBFWidth,etaBins,WMS[iList][iCanvas],WMSE[iList][iCanvas]); 
503               GetIntegral(gbfs[iList][iCanvas][a],integS[iList][iCanvas],integSE[iList][iCanvas]); 
504             }
505             else if(a==6){
506               GetWeightedMean(gbfs[iList][iCanvas][a],fStartBinBFWidthPhi,-1,WMPS[iList][iCanvas],WMPSE[iList][iCanvas]); // for phi calculate width 
507               GetIntegral(gbfs[iList][iCanvas][a],integPS[iList][iCanvas],integPSE[iList][iCanvas]); 
508             }
509           }
510           else{
511             fHistPNS[iList][a]->SetLineColor(2);
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();
515             fHistPPS[iList][a]->SetLineColor(1);
516             fHistPPS[iList][a]->ProjectionY(Form("pps%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
517             fHistNPS[iList][a]->SetLineColor(4);
518             fHistNPS[iList][a]->ProjectionY(Form("nps%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
519             fHistNNS[iList][a]->SetLineColor(8);
520             fHistNNS[iList][a]->ProjectionY(Form("nns%d",a),centralityArray[iCanvas]+1,centralityArray[iCanvas+1])->DrawCopy("same");
521           }
522         }
523       }
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       }
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]);
565       if(iList==0) gWM[i]->SetName("gCentrality");  
566       else gWM[i]->SetName(Form("gCentrality_%d",i));
567       gWMS[i] = new TGraphErrors(iCanvas,cent,WMS[i],centE,WMSE[i]); 
568       if(iList==0) gWMS[i]->SetName("gCentralityS");  
569       else gWMS[i]->SetName(Form("gCentralityS_%d",i)); 
570
571       gWMP[i] = new TGraphErrors(iCanvas,cent,WMP[i],centE,WMPE[i]);
572       if(iList==0) gWMP[i]->SetName("gCentralityPhi");  
573       else gWMP[i]->SetName(Form("gCentralityPhi_%d",i));
574       gWMPS[i] = new TGraphErrors(iCanvas,cent,WMPS[i],centE,WMPSE[i]); 
575       if(iList==0) gWMPS[i]->SetName("gCentralityPhiS");  
576       else gWMPS[i]->SetName(Form("gCentralityPhiS_%d",i)); 
577
578       ginteg[i] = new TGraphErrors(iCanvas,cent,integ[i],centE,integE[i]);
579       if(iList==0) ginteg[i]->SetName("gIntegral");  
580       else ginteg[i]->SetName(Form("gIntegral_%d",i));
581       gintegS[i] = new TGraphErrors(iCanvas,cent,integS[i],centE,integSE[i]); 
582       if(iList==0) gintegS[i]->SetName("gIntegralS");  
583       else gintegS[i]->SetName(Form("gIntegralS_%d",i)); 
584
585       gintegP[i] = new TGraphErrors(iCanvas,cent,integP[i],centE,integPE[i]);
586       if(iList==0) gintegP[i]->SetName("gIntegraPhil");  
587       else gintegP[i]->SetName(Form("gIntegralPhi_%d",i));
588       gintegPS[i] = new TGraphErrors(iCanvas,cent,integPS[i],centE,integPSE[i]); 
589       if(iList==0) gintegPS[i]->SetName("gIntegralPhiS");  
590       else gintegPS[i]->SetName(Form("gIntegralPhiS_%d",i)); 
591     }
592   }
593
594   TFile *fOut = NULL;
595   if(etaWindow > 0 || etaBins > -1){
596     if(correctWithEfficiency){
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       }
603     }
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     }
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     }
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   }
614   fOut->cd();
615   for(Int_t i = 0; i < iList+1; i++){
616     if(gbf[i][0][1]){
617     cout<<"PROCESS LIST "<<i<<" NOW!"<<endl;
618     for(Int_t a = 0; a < 7; a++){
619       if(a==1||a==6){
620       cout<<"PROCESS VARIABLE "<<a<<" NOW!"<<endl;
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       }
639
640       //printout in text format for delta eta      
641       for(Int_t j = 0; j < iCanvas; j++){
642         cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
643         
644         if(gbf[i][j][a]){
645           if(a==1){
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;
652           } 
653           else if(a==6){
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;
660           } 
661           gbf[i][j][a]->Write();
662           gbf[i][j][a]->Delete();
663         }
664         if(gbfs[i][j][a]){
665           if(a==1){
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;
672           } 
673           else if(a==6){
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;
680           } 
681           gbfs[i][j][a]->Write();
682           gbfs[i][j][a]->Delete();
683         }
684         cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
685         cout<<endl;
686       }
687       }
688     }
689
690     Double_t x,y;
691
692     cout<<"//================================ALICE================================//"<<endl;
693     if(gWM[i]){
694       cout<<"Double_t gWeightedMeanInEtaAlice[nCentralityBins] = {";
695       for(Int_t k = 0; k < gWM[i]->GetN()-1;k++){
696         gWM[i]->GetPoint(k,x,y);
697         cout<<y<<", ";      
698       }
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++){
704         cout<<gWM[i]->GetErrorY(k)<<", ";      
705       }
706       cout<<gWM[i]->GetErrorY(k)<<"};"<<endl;
707       gWM[i]->Write();
708     } 
709     if(gWMS[i]){
710       cout<<"Double_t gShuffledWeightedMeanInEtaAlice[nCentralityBins] = {";
711       for(Int_t k = 0; k < gWMS[i]->GetN()-1;k++){
712         gWMS[i]->GetPoint(k,x,y);
713         cout<<y<<", ";      
714       }
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++){
720         cout<<gWMS[i]->GetErrorY(k)<<", ";      
721       }
722       cout<<gWMS[i]->GetErrorY(k)<<"};"<<endl;
723       cout<<endl;
724       gWMS[i]->Write();
725     } 
726
727     if(gWMP[i]){
728       cout<<"Double_t gWeightedMeanInPhiAlice[nCentralityBins] = {";
729       for(Int_t k = 0; k < gWMP[i]->GetN()-1;k++){
730         gWMP[i]->GetPoint(k,x,y);
731         cout<<y<<", ";      
732       }
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++){
738         cout<<gWMP[i]->GetErrorY(k)<<", ";      
739       }
740       cout<<gWMP[i]->GetErrorY(k)<<"};"<<endl;
741       gWMP[i]->Write();
742     } 
743     if(gWMPS[i]){
744       cout<<"Double_t gShuffledWeightedMeanInPhiAlice[nCentralityBins] = {";
745       for(Int_t k = 0; k < gWMPS[i]->GetN()-1;k++){
746         gWMPS[i]->GetPoint(k,x,y);
747         cout<<y<<", ";      
748       }
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++){
754         cout<<gWMPS[i]->GetErrorY(k)<<", ";      
755       }
756       cout<<gWMPS[i]->GetErrorY(k)<<"};"<<endl;
757       cout<<endl;
758       gWMPS[i]->Write();
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();
796
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++){
817         gintegPS[i]->GetPoint(k,x,y);
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     }
834     }
835   }
836   fOut->Close();
837   f->Close();
838   gROOT->Reset();
839 }
840
841 //____________________________________________________________________//
842 void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME) {
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();
853   if(fStopBin > -1) fNumberOfBins = fStopBin;
854   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
855   
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;
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   
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;
887
888   WM  = delta;
889   WME = deltaError;
890 }
891
892 //____________________________________________________________________//
893 void 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   
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;
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   
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;
940
941   integ  = integral;
942   integE = integralError;
943 }
944
945 //___________________________________________________________//
946 // NOT USED any more
947 //___________________________________________________________//
948 void 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 }