8f78487b95bdbe31a58f3db9a732f51046abecd3
[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,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,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, 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 
98   Double_t WMPE[13][10];    // error
99   Double_t WMPS[13][10];     // weighted mean for phi (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_")){
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(5);
262           fHistN[iList][a]->RebinY(5);
263           fHistPP[iList][a]->RebinY(5);
264           fHistPN[iList][a]->RebinY(5);
265           fHistNP[iList][a]->RebinY(5);
266           fHistNN[iList][a]->RebinY(5);
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 (from 0.1 only!)
304               GetIntegral(gbf[iList][iCanvas][a],integ[iList][iCanvas],integE[iList][iCanvas]);
305             }
306             else if(a==6){
307               GetWeightedMean(gbf[iList][iCanvas][a],1,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(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(5);
352           fHistNS[iList][a]->RebinY(5);
353           fHistPPS[iList][a]->RebinY(5);
354           fHistPNS[iList][a]->RebinY(5);
355           fHistNPS[iList][a]->RebinY(5);
356           fHistNNS[iList][a]->RebinY(5);
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],1,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_%s", fStartBinBFWidth, fRebin,inFile.Data()),"RECREATE");
472   fOut->cd();
473   for(Int_t i = 0; i < iList+1; i++){
474     for(Int_t a = 0; a < 7; a++){
475       
476       if(fHistPN[i][a]){
477         (fHistPN[i][a]->ProjectionY(Form("hPN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
478         (fHistPP[i][a]->ProjectionY(Form("hPP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
479         (fHistNP[i][a]->ProjectionY(Form("hNP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
480         (fHistNN[i][a]->ProjectionY(Form("hNN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
481         (fHistP[i][a]->ProjectionY(Form("hP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
482         (fHistN[i][a]->ProjectionY(Form("hN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
483       }
484
485       if(fHistPNS[i][a]){
486         (fHistPNS[i][a]->ProjectionY(Form("hPNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
487         (fHistPPS[i][a]->ProjectionY(Form("hPPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
488         (fHistNPS[i][a]->ProjectionY(Form("hNPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
489         (fHistNNS[i][a]->ProjectionY(Form("hNNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
490         (fHistPS[i][a]->ProjectionY(Form("hPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
491         (fHistNS[i][a]->ProjectionY(Form("hNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
492       }
493       
494       for(Int_t j = 0; j < iCanvas; j++){
495         
496         if(gbf[i][j][a]){
497           //printout in text format for delta eta
498           if(a==1){
499             cout<<"Balance Function "<<i<<" "<<j<<" : ";
500             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
501             cout<<endl;
502             cout<<"Balance Function error "<<i<<" "<<j<<" : ";
503             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
504             cout<<endl;
505           } 
506           else if(a==6){
507             cout<<"Balance Function Phi "<<i<<" "<<j<<" : ";
508             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
509             cout<<endl;
510             cout<<"Balance Function Phi error "<<i<<" "<<j<<" : ";
511             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
512             cout<<endl;
513           } 
514           gbf[i][j][a]->Write();
515           gbf[i][j][a]->Delete();
516         }
517         if(gbfs[i][j][a]){
518           if(a==1){
519             cout<<"Balance Function (shuffled) "<<i<<" "<<j<<" : ";
520             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
521             cout<<endl;
522             cout<<"Balance Function error (shuffled) "<<i<<" "<<j<<" : ";
523             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
524             cout<<endl;
525           } 
526           else if(a==6){
527             cout<<"Balance Function Phi (shuffled) "<<i<<" "<<j<<" : ";
528             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
529             cout<<endl;
530             cout<<"Balance Function Phi error (shuffled) "<<i<<" "<<j<<" : ";
531             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
532             cout<<endl;
533           } 
534           gbfs[i][j][a]->Write();
535           gbfs[i][j][a]->Delete();
536         }
537       }
538     }
539
540     Double_t x,y;
541       
542     if(gWM[i]){
543       cout<<"Balance Function WM "<<i<<" : ";
544       for(Int_t k = 0; k < gWM[i]->GetN();k++){
545         gWM[i]->GetPoint(k,x,y);
546         cout<<y<<", ";      
547       }
548       cout<<endl;
549       cout<<"Balance Function WM error "<<i<<" : ";
550       for(Int_t k = 0; k < gWM[i]->GetN();k++){
551         cout<<gWM[i]->GetErrorY(k)<<", ";      
552       }
553       cout<<endl;
554       gWM[i]->Write();
555     } 
556     if(gWMS[i]){
557       cout<<"Balance Function WM (shuffled) "<<i<<" : ";
558       for(Int_t k = 0; k < gWMS[i]->GetN();k++){
559         gWMS[i]->GetPoint(k,x,y);
560         cout<<y<<", ";      
561       }
562       cout<<endl;
563       cout<<"Balance Function WM error (shuffled) "<<i<<" : ";
564       for(Int_t k = 0; k < gWMS[i]->GetN();k++){
565         cout<<gWMS[i]->GetErrorY(k)<<", ";      
566       }
567       cout<<endl;
568       gWMS[i]->Write();
569     } 
570
571    if(gWMP[i]){
572       cout<<"Balance Function Phi WMP "<<i<<" : ";
573       for(Int_t k = 0; k < gWMP[i]->GetN();k++){
574         gWMP[i]->GetPoint(k,x,y);
575         cout<<y<<", ";      
576       }
577       cout<<endl;
578       cout<<"Balance Function Phi WMP error "<<i<<" : ";
579       for(Int_t k = 0; k < gWMP[i]->GetN();k++){
580         cout<<gWMP[i]->GetErrorY(k)<<", ";      
581       }
582       cout<<endl;
583       gWMP[i]->Write();
584     } 
585     if(gWMPS[i]){
586       cout<<"Balance Function Phi WMP (shuffled) "<<i<<" : ";
587       for(Int_t k = 0; k < gWMPS[i]->GetN();k++){
588         gWMPS[i]->GetPoint(k,x,y);
589         cout<<y<<", ";      
590       }
591       cout<<endl;
592       cout<<"Balance Function Phi WMP error (shuffled) "<<i<<" : ";
593       for(Int_t k = 0; k < gWMPS[i]->GetN();k++){
594         cout<<gWMPS[i]->GetErrorY(k)<<", ";      
595       }
596       cout<<endl;
597       gWMPS[i]->Write();
598     } 
599
600     if(ginteg[i]) ginteg[i]->Write();
601     if(gintegS[i]) gintegS[i]->Write();
602     if(gintegP[i]) gintegP[i]->Write();
603     if(gintegPS[i]) gintegPS[i]->Write();
604   }
605   fOut->Close();
606   f->Close();
607   gROOT->Reset();
608 }
609
610 //____________________________________________________________________//
611 void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Double_t &WM, Double_t &WME) {
612
613   //Prints the calculated width of the BF and its error
614   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
615   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
616   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
617   Double_t deltaBalP2 = 0.0, integral = 0.0;
618   Double_t deltaErrorNew = 0.0;
619
620   //Retrieve this variables from Histogram
621   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
622   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
623   
624   cout<<"=================================================="<<endl;
625   cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
626   cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
627   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
628     cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
629   } 
630   cout<<"=================================================="<<endl;
631   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
632     gSumXi += gHistBalance->GetBinCenter(i);
633     gSumBi += gHistBalance->GetBinContent(i);
634     gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
635     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
636     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
637     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
638     gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
639     
640     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
641     integral += fP2Step*gHistBalance->GetBinContent(i);
642   }
643   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
644     deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
645   
646   Double_t integralError = TMath::Sqrt(deltaBalP2);
647   
648   Double_t delta = gSumBiXi / gSumBi;
649   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
650   
651   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
652   cout<<"New error: "<<deltaErrorNew<<endl;
653   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
654   cout<<"=================================================="<<endl;
655
656   WM  = delta;
657   WME = deltaError;
658 }
659
660 //____________________________________________________________________//
661 void GetIntegral(TH1D *gHistBalance, Double_t &integ, Double_t &integE) {
662
663   //always start with 1
664   Int_t fStartBin = 1; 
665
666   //Prints the calculated width of the BF and its error
667   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
668   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
669   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
670   Double_t deltaBalP2 = 0.0, integral = 0.0;
671   Double_t deltaErrorNew = 0.0;
672
673   //Retrieve this variables from Histogram
674   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
675   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
676   
677   cout<<"=================================================="<<endl;
678   cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
679   cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
680   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
681     cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
682   } 
683   cout<<"=================================================="<<endl;
684   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
685     gSumXi += gHistBalance->GetBinCenter(i);
686     gSumBi += gHistBalance->GetBinContent(i);
687     gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
688     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
689     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
690     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
691     gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
692     
693     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
694     integral += fP2Step*gHistBalance->GetBinContent(i);
695   }
696   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
697     deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
698   
699   Double_t integralError = TMath::Sqrt(deltaBalP2);
700   
701   Double_t delta = gSumBiXi / gSumBi;
702   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
703   
704   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
705   cout<<"New error: "<<deltaErrorNew<<endl;
706   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
707   cout<<"=================================================="<<endl;
708
709   integ  = integral;
710   integE = integralError;
711 }
712
713 //___________________________________________________________//
714 // NOT USED any more
715 //___________________________________________________________//
716 void mergeOutput(const char* outputDir) {
717   //Function to merge the output of the sub-jobs
718   //Create a BF object
719   AliBalance *bf = new AliBalance();
720
721   //connect to AliEn's API services
722   TGrid::Connect("alien://"); 
723
724   //Getting the output dir from the env. variable 
725   //(JDL field: JDLVariables={"OutputDir"};)
726   TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");
727   
728   Int_t nEntries = result->GetEntries();
729
730   TString alienUrl;
731   TDirectoryFile *dirSubJob;
732
733   TString gCutName[4] = {"Total","Offline trigger",
734                          "Vertex","Analyzed"};
735   TH1F *fHistEventStats = new TH1F("fHistEventStats",
736                                    "Event statistics;;N_{events}",
737                                    4,0.5,4.5);
738   for(Int_t i = 1; i <= 4; i++)
739     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
740
741   AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");
742   for(Int_t i = 0; i < nEntries; i++) {
743     alienUrl = result->GetKey(i,"turl");
744     alienUrl += "#AnalysisResults.root";
745     Printf("Opening file: %s",alienUrl.Data());
746     TFile *file = TFile::Open(alienUrl.Data());
747     dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWGCFEbyE.outputBalanceFunctionAnalysis.root"));
748
749     //merge BF
750     AliBalance *bfSubJob = dynamic_cast<AliBalance *>(dirSubJob->Get("AliBalance"));
751     //bfSubJob->PrintResults();
752     bf->Merge(bfSubJob);
753     //delete bfSubJob;
754
755     //merge event stats
756     TList *listSubJob = dynamic_cast<TList *>(dirSubJob->Get("listQA"));
757     fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));
758
759     bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));
760     delete listSubJob;
761   }
762
763   //Create the output file
764   TString outputFile = "AnalysisResults.Merged.root";
765   TFile *foutput = TFile::Open(outputFile.Data(),"recreate");
766   TDirectoryFile *dirOutput = new TDirectoryFile();
767   dirOutput->SetName("PWGCFEbyE.outputBalanceFunctionAnalysis.root");
768   //dirOutput->cd();
769   dirOutput->Add(bf);
770   TList *list = new TList();
771   list->SetName("listQA");
772   list->Add(fHistEventStats);
773   list->Add(bfTrackCuts);
774   dirOutput->Add(list);
775   dirOutput->Write();
776   bf->Write();
777   list->Write();
778   foutput->Close();
779
780     //cout<<alienUrl.Data()<<endl;
781 }