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