Added optional Acceptance correction for Balance Functions in Delta Eta
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 24 Jan 2012 21:07:41 +0000 (21:07 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 24 Jan 2012 21:07:41 +0000 (21:07 +0000)
PWGCF/EBYE/BalanceFunctions/AliBalance.cxx
PWGCF/EBYE/BalanceFunctions/AliBalance.h
PWGCF/EBYE/macros/readBalanceFunction.C

index 431d8d0..30b6575 100644 (file)
@@ -701,9 +701,15 @@ void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
 }
  
 //____________________________________________________________________//
-TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax) {
+TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow) {
   //Returns the BF histogram, extracted from the 6 TH2D objects 
   //(private members) of the AliBalance class.
+  //
+  // Acceptance correction: 
+  // - only for analysis type = kEta
+  // - only if etaWindow > 0 (default = -1.)
+  // - calculated as proposed by STAR 
+  //
   TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
   TString histName = "gHistBalanceFunctionHistogram";
   histName += gAnalysisType[iAnalysisType];
@@ -774,6 +780,17 @@ TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centr
     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
     gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
   }
+
+  // do the acceptance correction (only for Eta and etaWindow > 0)
+  if(iAnalysisType == kEta && etaWindow > 0){
+    for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
+      
+      Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
+      Double_t corrected    = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow );
+      gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected);
+      
+    }
+  }
   
   PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
 
index aaddc20..679ea53 100644 (file)
@@ -106,7 +106,7 @@ class AliBalance : public TObject {
   void SetHistNnn(Int_t iAnalysisType, TH2D *gHist) { 
     fHistNN[iAnalysisType] = gHist;}
 
-  TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax);
+  TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow = -1);
   void PrintResults(Int_t iAnalysisType, TH1D *gHist);
 
  private:
index 6e661eb..8f78487 100644 (file)
-const TString gBFAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};\r
-\r
-const Int_t nrOfCentralities = 9;\r
-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)\r
-//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)\r
-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 \r
-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)\r
-\r
-void readBalanceFunction(Bool_t bHistos = kTRUE, TString inFile = "AnalysisResults.root",Int_t fStartBinBFWidth = 1, Int_t fRebin = 1,TString centEst = "V0M",TString extraString = "") {\r
-  // Macro to read the output of the BF analysis:  MW: CHANGE THIS!!!!\r
-  //i) Prints and draws the final BF output\r
-  //ii) Plots the QA part of the analysis\r
-  //iii) store BF in output file\r
-  //Author: Panos.Christakoglou@cern.ch, m.weber@cern.ch\r
-  //Loading the needed libraries\r
-  gSystem->Load("libProofPlayer.so");\r
-  gSystem->Load("libANALYSIS.so");\r
-  gSystem->Load("libANALYSISalice.so");\r
-  gSystem->Load("libPWGCFebye.so");\r
-\r
-  //Draw BF              \r
-  drawBF(bHistos,inFile, fStartBinBFWidth, fRebin,centEst,extraString);\r
-\r
-  //Merge the output\r
-  //mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");\r
-}\r
-\r
-//___________________________________________________________//\r
-void drawBF(Bool_t bHistos = kTRUE, TString inFile = "AnalysisResults.root", Int_t fStartBinBFWidth = 1, Int_t fRebin = 1, TString centEst = "V0M",TString extraString = "") {\r
-  //Function to draw the BF objects and write them into the output file\r
-\r
-  Int_t maximumCanvases = 10;\r
-  Int_t iCanvas         = 0;\r
-  TCanvas *cQA[10];\r
-  TCanvas *cQAV0M = new TCanvas("cQAV0M","V0M multiplicities");\r
-  cQAV0M->Divide(4,3);\r
-  TCanvas *cQARef = new TCanvas("cQARef","reference track multiplicities");\r
-  cQARef->Divide(4,3);\r
-  TCanvas *cBF[10];\r
-  TCanvas *cBFS[10];\r
-\r
-  // get the file\r
-  TFile *f = TFile::Open(inFile.Data());\r
-  if(!f) {\r
-    Printf("File not found!!!");\r
-    break;\r
-  }\r
-\r
-  // get the BF output directory\r
-  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis"));\r
-  if(!dir) {\r
-    Printf("Output directory not found!!!");\r
-    break;\r
-  }\r
-\r
-  // loop over all lists and plot the BF and QA\r
-  TList *list = NULL;\r
-  TString listName;\r
-  TIter nextkey( dir->GetListOfKeys() );\r
-  TKey *key;\r
-\r
-  AliBalance *bf[7];\r
-  AliBalance *bfs[7];\r
-  TH1D *gbf[10][7];\r
-  TH1D *gbfs[10][7];\r
-\r
-  for(Int_t i = 0; i < 10; i++){\r
-    for(Int_t j = 0; j < 7; j++){\r
-      gbf[i][j] = NULL;\r
-      gbfs[i][j] = NULL;\r
-    }\r
-  }\r
-\r
-  TH2D *fHistP[7]; //N+\r
-  TH2D *fHistN[7]; //N-\r
-  TH2D *fHistPN[7]; //N+-\r
-  TH2D *fHistNP[7]; //N-+\r
-  TH2D *fHistPP[7]; //N++\r
-  TH2D *fHistNN[7]; //N--\r
-\r
-  TH2D *fHistPS[7]; //N+\r
-  TH2D *fHistNS[7]; //N-\r
-  TH2D *fHistPNS[7]; //N+-\r
-  TH2D *fHistNPS[7]; //N-+\r
-  TH2D *fHistPPS[7]; //N++\r
-  TH2D *fHistNNS[7]; //N--\r
-\r
-  Double_t WM[10];     // weighted mean for eta (recalculated from fStartBin)\r
-  Double_t WME[10];    // error\r
-  Double_t WMS[10];     // weighted mean for eta (recalculated from fStartBin) (shuffled)\r
-  Double_t WMSE[10];    // error (shuffled)\r
-\r
-  for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){\r
-    \r
-    cBF[iCanvas] = new TCanvas(Form("cBF_%d",iCanvas),Form("cBF_%d",iCanvas));\r
-    cBF[iCanvas]->Divide(3,3);\r
-    \r
-    cBFS[iCanvas] = new TCanvas(Form("Shuffled_%d",iCanvas),Form("Shuffled_%d",iCanvas));\r
-    cBFS[iCanvas]->Divide(3,3);\r
-    \r
-  }\r
-\r
-  while ( (key = (TKey*)nextkey())) {\r
-\r
-    list = (TList*)key->ReadObj();\r
-    listName = TString(list->GetName());\r
-\r
-    cout<<"Processing list "<<listName<<endl;\r
-\r
-    // ----------------------------------------------------\r
-    // plot QA histograms\r
-    if(listName.Contains("QA")&&listName.Contains(extraString.Data())){     \r
-\r
-      cQA[iCanvas] = new TCanvas(listName,listName);\r
-      cQA[iCanvas]->Divide(4,3);\r
-\r
-     cQA[iCanvas]->cd(1);\r
-      TH1F* histEventStats = (TH1F*)list->FindObject("fHistEventStats");\r
-      if(histEventStats){\r
-       histEventStats->SetFillColor(9);\r
-       histEventStats->Draw();\r
-      }\r
-\r
-      cQA[iCanvas]->cd(2);\r
-      TH1F* histTriggerStats = (TH1F*)list->FindObject("fHistTriggerStats");\r
-      if(histTriggerStats){\r
-       histTriggerStats->SetFillColor(9);\r
-       histTriggerStats->Draw();\r
-      }\r
-\r
-     cQA[iCanvas]->cd(3);\r
-      TH1F* histTrackStats = (TH1F*)list->FindObject("fHistTrackStats");\r
-      if(histTrackStats){\r
-       histTrackStats->SetFillColor(9);\r
-       histTrackStats->Draw();\r
-      }\r
-\r
-      cQA[iCanvas]->cd(4);\r
-      TH1F* histCentStats = (TH1F*)list->FindObject("fHistCentStats");\r
-      if(histCentStats){\r
-       histCentStats->SetFillColor(9);\r
-       histCentStats->Draw("colz");\r
-      }\r
-\r
-\r
-\r
-      cQA[iCanvas]->cd(5);\r
-      TH1F* histVx = (TH1F*)list->FindObject("fHistVx");\r
-      if(histVx){\r
-       histVx->SetFillColor(9);\r
-       histVx->Draw();\r
-      }\r
-      cQA[iCanvas]->cd(6);\r
-      TH1F* histVy = (TH1F*)list->FindObject("fHistVy");\r
-      if(histVy){\r
-       histVy->SetFillColor(9);\r
-       histVy->Draw();\r
-      }\r
-      \r
-      cQA[iCanvas]->cd(7);\r
-      TH1F* histVz = (TH1F*)list->FindObject("fHistVz");\r
-      if(histVz){\r
-       histVz->SetFillColor(9);\r
-       histVz->Draw();\r
-      }\r
-\r
-      cQA[iCanvas]->cd(8);\r
-      cQA[iCanvas]->cd(8)->SetLogz();\r
-      TH2F* histDCA = (TH2F*)list->FindObject("fHistDCA");  \r
-      if(histDCA) histDCA->Draw("colz");\r
-      \r
-\r
-      cQA[iCanvas]->cd(9);\r
-      cQA[iCanvas]->cd(9)->SetLogz();\r
-      TH2F* histClus = (TH2F*)list->FindObject("fHistClus");\r
-      if(histClus) histClus->Draw("colz");\r
-      \r
-      cQA[iCanvas]->cd(10);\r
-      TH1F* histPt = (TH1F*)list->FindObject("fHistPt");\r
-      if(histPt){\r
-       histPt->SetFillColor(9);\r
-       histPt->Draw();\r
-      }\r
-      \r
-      cQA[iCanvas]->cd(11);\r
-      TH1F* histEta = (TH1F*)list->FindObject("fHistEta");\r
-      if(histEta){\r
-       histEta->SetFillColor(9);\r
-       histEta->Draw();\r
-      }\r
-      \r
-      cQA[iCanvas]->cd(12);\r
-      TH1F* histPhi = (TH1F*)list->FindObject("fHistPhi");\r
-      if(histPhi){\r
-       histPhi->SetFillColor(9);\r
-       histPhi->Draw();\r
-      }\r
-\r
-      // centrality estimator QA\r
-      cQAV0M->cd(iCanvas+1);\r
-      cQAV0M->cd(iCanvas+1)->SetLogz();\r
-      TH1F* histV0M = (TH1F*)list->FindObject("fHistV0M");\r
-      if(histV0M){\r
-       histV0M->Draw("colz");\r
-      }\r
-\r
-      cQARef->cd(iCanvas+1);\r
-      cQARef->cd(iCanvas+1)->SetLogz();\r
-      TH1F* histRef = (TH1F*)list->FindObject("fHistRefTracks");\r
-      if(histRef){\r
-       histRef->Draw("colz");\r
-      }\r
-    }\r
-    // ----------------------------------------------------\r
-\r
-    // ----------------------------------------------------\r
-    // calculate and plot BF \r
-    if(listName.Contains("BF_")&&listName.Contains(extraString.Data())){\r
-\r
-\r
-      for(Int_t a = 0; a < 7; a++){\r
-\r
-       cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;\r
-\r
-       // create the BF object\r
-       bf[a]  = new AliBalance();\r
-\r
-       fHistP[a] = (TH2D*)list->FindObject(Form("fHistP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistN[a] = (TH2D*)list->FindObject(Form("fHistN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistPP[a] = (TH2D*)list->FindObject(Form("fHistPP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistPN[a] = (TH2D*)list->FindObject(Form("fHistPN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistNP[a] = (TH2D*)list->FindObject(Form("fHistNP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistNN[a] = (TH2D*)list->FindObject(Form("fHistNN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-\r
-       // rebin histograms (be careful with divider!)\r
-       fHistP[a]->RebinY(fRebin);\r
-       fHistN[a]->RebinY(fRebin);\r
-       fHistPP[a]->RebinY(fRebin);\r
-       fHistPN[a]->RebinY(fRebin);\r
-       fHistNP[a]->RebinY(fRebin);\r
-       fHistNN[a]->RebinY(fRebin);\r
-\r
-       fHistP[a]->SetName(Form("%s_%d",fHistP[a]->GetName(),iCanvas));\r
-       fHistN[a]->SetName(Form("%s_%d",fHistN[a]->GetName(),iCanvas));\r
-       fHistPP[a]->SetName(Form("%s_%d",fHistPP[a]->GetName(),iCanvas));\r
-       fHistPN[a]->SetName(Form("%s_%d",fHistPN[a]->GetName(),iCanvas));\r
-       fHistNP[a]->SetName(Form("%s_%d",fHistNP[a]->GetName(),iCanvas));\r
-       fHistNN[a]->SetName(Form("%s_%d",fHistNN[a]->GetName(),iCanvas));\r
-\r
-       // set histograms in AliBalance object\r
-       bf[a]->SetHistNp(a, fHistP[a]);\r
-       bf[a]->SetHistNn(a, fHistN[a]);\r
-       bf[a]->SetHistNpp(a, fHistPP[a]);\r
-       bf[a]->SetHistNpn(a, fHistPN[a]);\r
-       bf[a]->SetHistNnp(a, fHistNP[a]);\r
-       bf[a]->SetHistNnn(a, fHistNN[a]);\r
-\r
-       for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){\r
-\r
-         gbf[iCanvas][a] = bf[a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1]);\r
-         gbf[iCanvas][a]->SetName(Form("BF_%s_Cent_%.0f_%.0f",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1]));\r
-\r
-         cBF[iCanvas]->cd(a+1);\r
-         gbf[iCanvas][a]->SetMarkerStyle(20);\r
-\r
-         if(!bHistos){\r
-           gbf[iCanvas][a]->DrawCopy("AP");\r
-           if(a==1) GetWeightedMean(gbf[iCanvas][a],fStartBinBFWidth,WM[iCanvas],WME[iCanvas]); // for eta recalculate width (from 0.1 only!)\r
-         }\r
-         else{\r
-           fHistPN[a]->SetLineColor(2);\r
-           fHistPN[a]->ProjectionY(Form("pn%d",a))->DrawCopy();\r
-           fHistPP[a]->SetLineColor(1);\r
-           fHistPP[a]->ProjectionY(Form("pp%d",a))->DrawCopy("same");\r
-           fHistNP[a]->SetLineColor(4);\r
-           fHistNP[a]->ProjectionY(Form("np%d",a))->DrawCopy("same");\r
-           fHistNN[a]->SetLineColor(8);\r
-           fHistNN[a]->ProjectionY(Form("nn%d",a))->DrawCopy("same");\r
-         }\r
-       }\r
-      }\r
-    }\r
-    // ----------------------------------------------------\r
-\r
-    // ----------------------------------------------------\r
-    // calculate and plot BF (shuffled)\r
-    if(listName.Contains("BFShuffled")&&listName.Contains(extraString.Data())){\r
-\r
-      for(Int_t a = 0; a < 7; a++){\r
-\r
-       // create the BF object\r
-       bfs[a]  = new AliBalance();\r
-\r
-       fHistPS[a] = (TH2D*)list->FindObject(Form("fHistP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistNS[a] = (TH2D*)list->FindObject(Form("fHistN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistPPS[a] = (TH2D*)list->FindObject(Form("fHistPP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistPNS[a] = (TH2D*)list->FindObject(Form("fHistPN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistNPS[a] = (TH2D*)list->FindObject(Form("fHistNP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-       fHistNNS[a] = (TH2D*)list->FindObject(Form("fHistNN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));\r
-\r
-       // rebin histograms (be careful with divider!)\r
-       fHistPS[a]->RebinY(fRebin);\r
-       fHistNS[a]->RebinY(fRebin);\r
-       fHistPPS[a]->RebinY(fRebin);\r
-       fHistPNS[a]->RebinY(fRebin);\r
-       fHistNPS[a]->RebinY(fRebin);\r
-       fHistNNS[a]->RebinY(fRebin);\r
-\r
-       fHistPS[a]->SetName(Form("%s_%d",fHistPS[a]->GetName(),iCanvas));\r
-       fHistNS[a]->SetName(Form("%s_%d",fHistNS[a]->GetName(),iCanvas));\r
-       fHistPPS[a]->SetName(Form("%s_%d",fHistPPS[a]->GetName(),iCanvas));\r
-       fHistPNS[a]->SetName(Form("%s_%d",fHistPNS[a]->GetName(),iCanvas));\r
-       fHistNPS[a]->SetName(Form("%s_%d",fHistNPS[a]->GetName(),iCanvas));\r
-       fHistNNS[a]->SetName(Form("%s_%d",fHistNNS[a]->GetName(),iCanvas));\r
-\r
-       // set histograms in AliBalance object\r
-       bfs[a]->SetHistNp(a, fHistPS[a]);\r
-       bfs[a]->SetHistNn(a, fHistNS[a]);\r
-       bfs[a]->SetHistNpp(a, fHistPPS[a]);\r
-       bfs[a]->SetHistNpn(a, fHistPNS[a]);\r
-       bfs[a]->SetHistNnp(a, fHistNPS[a]);\r
-       bfs[a]->SetHistNnn(a, fHistNNS[a]);\r
-\r
-       for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){\r
-         \r
-         gbfs[iCanvas][a] = bfs[a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1]);\r
-         gbf[iCanvas][a]->SetName(Form("BFS_%s_Cent_%.0f_%.0f",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1]));\r
-         \r
-         cBFS[iCanvas]->cd(a+1);\r
-         gbfs[iCanvas][a]->SetMarkerStyle(20);\r
-         if(!bHistos){\r
-           gbfs[iCanvas][a]->DrawCopy("AP");\r
-           if(a==1) GetWeightedMean(gbfs[iCanvas][a],fStartBinBFWidth,WMS[iCanvas],WMSE[iCanvas]); // for eta recalculate width (from 0.1 only!)\r
-         }\r
-         else{\r
-           fHistPNS[a]->SetLineColor(2);\r
-           fHistPNS[a]->ProjectionY(Form("pns%d",a))->DrawCopy();\r
-           fHistPPS[a]->SetLineColor(1);\r
-           fHistPPS[a]->ProjectionY(Form("pps%d",a))->DrawCopy("same");\r
-           fHistNPS[a]->SetLineColor(4);\r
-           fHistNPS[a]->ProjectionY(Form("nps%d",a))->DrawCopy("same");\r
-           fHistNNS[a]->SetLineColor(8);\r
-           fHistNNS[a]->ProjectionY(Form("nns%d",a))->DrawCopy("same");\r
-         }\r
-       }\r
-      }\r
-    }\r
-    // ----------------------------------------------------\r
-  }\r
-\r
-  // for BF calculation create also graphs with weighted mean for eta\r
-  TGraphErrors *gWM  = NULL;\r
-  TGraphErrors *gWMS = NULL;\r
-  if(!bHistos && gbf[1][1]){\r
-    gWM = new TGraphErrors(iCanvas,cent,WM,centE,WME);\r
-    gWM->SetName("gCentrality");    \r
-  }\r
-  if(!bHistos && gbfs[1][1]){\r
-    gWMS = new TGraphErrors(iCanvas,cent,WMS,centE,WMSE);  \r
-    gWMS->SetName("gCentralityS");\r
-  }\r
-\r
-  TFile *fOut = TFile::Open(Form("Histograms_WMstart%d_rebin%d_%s", fStartBinBFWidth, fRebin,inFile.Data()),"RECREATE");\r
-  fOut->cd();\r
-  for(Int_t a = 0; a < 7; a++){\r
-\r
-    if(fHistPN[a]){\r
-      (fHistPN[a]->ProjectionY(Form("hPN_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistPP[a]->ProjectionY(Form("hPP_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistNP[a]->ProjectionY(Form("hNP_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistNN[a]->ProjectionY(Form("hNN_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistP[a]->ProjectionY(Form("hP_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistN[a]->ProjectionY(Form("hN_%s",gBFAnalysisType[a].Data())))->Write();\r
-    }\r
-\r
-    if(fHistPNS[a]){\r
-      (fHistPNS[a]->ProjectionY(Form("hPNS_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistPPS[a]->ProjectionY(Form("hPPS_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistNPS[a]->ProjectionY(Form("hNPS_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistNNS[a]->ProjectionY(Form("hNNS_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistPS[a]->ProjectionY(Form("hPS_%s",gBFAnalysisType[a].Data())))->Write();\r
-      (fHistNS[a]->ProjectionY(Form("hNS_%s",gBFAnalysisType[a].Data())))->Write();\r
-    }\r
-\r
-    for(Int_t i = 0; i < iCanvas; i++){\r
-      \r
-      if(gbf[i][a]){\r
-       gbf[i][a]->Write();\r
-       gbf[i][a]->Delete();\r
-      }\r
-      if(gbfs[i][a]){\r
-       gbfs[i][a]->Write();\r
-       gbfs[i][a]->Delete();\r
-      }\r
-    }\r
-  }\r
-\r
-  if(gWM) gWM->Write();\r
-  if(gWMS) gWMS->Write();\r
-\r
-  fOut->Close();\r
-  f->Close();\r
-  gROOT->Reset();\r
-}\r
-\r
-//____________________________________________________________________//\r
-void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Double_t &WM, Double_t &WME) {\r
-\r
-  //Prints the calculated width of the BF and its error\r
-  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;\r
-  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;\r
-  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;\r
-  Double_t deltaBalP2 = 0.0, integral = 0.0;\r
-  Double_t deltaErrorNew = 0.0;\r
-\r
-  //Retrieve this variables from Histogram\r
-  Int_t fNumberOfBins = gHistBalance->GetNbinsX();\r
-  Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!\r
-  \r
-  cout<<"=================================================="<<endl;\r
-  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;\r
-  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;\r
-  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { \r
-    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;\r
-  } \r
-  cout<<"=================================================="<<endl;\r
-  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {\r
-    gSumXi += gHistBalance->GetBinCenter(i);\r
-    gSumBi += gHistBalance->GetBinContent(i);\r
-    gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);\r
-    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);\r
-    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);\r
-    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);\r
-    gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);\r
-    \r
-    deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);\r
-    integral += fP2Step*gHistBalance->GetBinContent(i);\r
-  }\r
-  for(Int_t i = 1; i < fNumberOfBins; i++)\r
-    deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);\r
-  \r
-  Double_t integralError = TMath::Sqrt(deltaBalP2);\r
-  \r
-  Double_t delta = gSumBiXi / gSumBi;\r
-  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );\r
-  \r
-  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;\r
-  cout<<"New error: "<<deltaErrorNew<<endl;\r
-  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;\r
-  cout<<"=================================================="<<endl;\r
-\r
-  WM  = delta;\r
-  WME = deltaError;\r
-}\r
-\r
-//___________________________________________________________//\r
-// NOT USED any more\r
-//___________________________________________________________//\r
-void mergeOutput(const char* outputDir) {\r
-  //Function to merge the output of the sub-jobs\r
-  //Create a BF object\r
-  AliBalance *bf = new AliBalance();\r
-\r
-  //connect to AliEn's API services\r
-  TGrid::Connect("alien://"); \r
-\r
-  //Getting the output dir from the env. variable \r
-  //(JDL field: JDLVariables={"OutputDir"};)\r
-  TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");\r
-  \r
-  Int_t nEntries = result->GetEntries();\r
-\r
-  TString alienUrl;\r
-  TDirectoryFile *dirSubJob;\r
-\r
-  TString gCutName[4] = {"Total","Offline trigger",\r
-                         "Vertex","Analyzed"};\r
-  TH1F *fHistEventStats = new TH1F("fHistEventStats",\r
-                                  "Event statistics;;N_{events}",\r
-                                  4,0.5,4.5);\r
-  for(Int_t i = 1; i <= 4; i++)\r
-    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
-\r
-  AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");\r
-  for(Int_t i = 0; i < nEntries; i++) {\r
-    alienUrl = result->GetKey(i,"turl");\r
-    alienUrl += "#AnalysisResults.root";\r
-    Printf("Opening file: %s",alienUrl.Data());\r
-    TFile *file = TFile::Open(alienUrl.Data());\r
-    dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWGCFEbyE.outputBalanceFunctionAnalysis.root"));\r
-\r
-    //merge BF\r
-    AliBalance *bfSubJob = dynamic_cast<AliBalance *>(dirSubJob->Get("AliBalance"));\r
-    //bfSubJob->PrintResults();\r
-    bf->Merge(bfSubJob);\r
-    //delete bfSubJob;\r
-\r
-    //merge event stats\r
-    TList *listSubJob = dynamic_cast<TList *>(dirSubJob->Get("listQA"));\r
-    fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));\r
-\r
-    bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));\r
-    delete listSubJob;\r
-  }\r
-\r
-  //Create the output file\r
-  TString outputFile = "AnalysisResults.Merged.root";\r
-  TFile *foutput = TFile::Open(outputFile.Data(),"recreate");\r
-  TDirectoryFile *dirOutput = new TDirectoryFile();\r
-  dirOutput->SetName("PWGCFEbyE.outputBalanceFunctionAnalysis.root");\r
-  //dirOutput->cd();\r
-  dirOutput->Add(bf);\r
-  TList *list = new TList();\r
-  list->SetName("listQA");\r
-  list->Add(fHistEventStats);\r
-  list->Add(bfTrackCuts);\r
-  dirOutput->Add(list);\r
-  dirOutput->Write();\r
-  bf->Write();\r
-  list->Write();\r
-  foutput->Close();\r
-\r
-    //cout<<alienUrl.Data()<<endl;\r
-}\r
+const TString gBFAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
+
+const Int_t nrOfCentralities = 9;
+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)
+//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)
+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 
+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)
+
+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) {
+  // Macro to read the output of the BF analysis:  MW: CHANGE THIS!!!!
+  //i) Prints and draws the final BF output
+  //ii) Plots the QA part of the analysis
+  //iii) store BF in output file
+  //Author: Panos.Christakoglou@cern.ch, m.weber@cern.ch
+  //Loading the needed libraries
+  gSystem->Load("libProofPlayer.so");
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libPWGCFebye.so");
+
+  //Draw BF       
+  drawBF(bHistos,inFile, fStartBinBFWidth, fRebin,centEst, "",  etaWindow);    
+  
+
+  //Merge the output
+  //mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");
+}
+
+//___________________________________________________________//
+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) {
+  //Function to draw the BF objects and write them into the output file
+
+  Int_t maximumCanvases = 13;
+  Int_t iCanvas         = 0;
+  Int_t iList           = -1;
+  TCanvas *cQA[13][10];
+  TCanvas *cQAV0M = new TCanvas("cQAV0M","V0M multiplicities");
+  cQAV0M->Divide(4,3);
+  TCanvas *cQARef = new TCanvas("cQARef","reference track multiplicities");
+  cQARef->Divide(4,3);
+  TCanvas *cBF[13][10];
+  TCanvas *cBFS[13][10];
+
+  // get the file
+  TFile *f = TFile::Open(inFile.Data());
+  if(!f) {
+    Printf("File not found!!!");
+    break;
+  }
+
+  // get the BF output directory
+  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionAnalysis"));
+  if(!dir) {
+    Printf("Output directory not found!!!");
+    break;
+  }
+
+  // loop over all lists and plot the BF and QA
+  TList *list = NULL;
+  TString listName;
+  TIter nextkey( dir->GetListOfKeys() );
+  TKey *key;
+
+  AliBalance *bf[13][7];
+  AliBalance *bfs[13][7];
+  TH1D *gbf[13][10][7];
+  TH1D *gbfs[13][10][7];
+
+  for(Int_t i = 0; i < 13; i++){
+    for(Int_t j = 0; j < 10; j++){
+      for(Int_t k = 0; k < 7; k++){
+       gbf[i][j][k]  = NULL;
+       gbfs[i][j][k] = NULL;
+      }
+    }
+  }
+
+  TH2D *fHistP[13][7]; //N+
+  TH2D *fHistN[13][7]; //N-
+  TH2D *fHistPN[13][7]; //N+-
+  TH2D *fHistNP[13][7]; //N-+
+  TH2D *fHistPP[13][7]; //N++
+  TH2D *fHistNN[13][7]; //N--
+
+  TH2D *fHistPS[13][7]; //N+
+  TH2D *fHistNS[13][7]; //N-
+  TH2D *fHistPNS[13][7]; //N+-
+  TH2D *fHistNPS[13][7]; //N-+
+  TH2D *fHistPPS[13][7]; //N++
+  TH2D *fHistNNS[13][7]; //N--
+
+  Double_t WM[13][10];     // weighted mean for eta (recalculated from fStartBin)
+  Double_t WME[13][10];    // error
+  Double_t WMS[13][10];     // weighted mean for eta (recalculated from fStartBin) (shuffled)
+  Double_t WMSE[13][10];    // error (shuffled)
+
+  Double_t WMP[13][10];     // weighted mean for phi 
+  Double_t WMPE[13][10];    // error
+  Double_t WMPS[13][10];     // weighted mean for phi (shuffled)
+  Double_t WMPSE[13][10];    // error (shuffled)
+
+  Double_t integ[13][10];     // integral for eta (calculated from bin 1)
+  Double_t integE[13][10];    // error
+  Double_t integS[13][10];     // integral for eta (calculated from bin 1) (shuffled)
+  Double_t integSE[13][10];    // error (shuffled)
+
+  Double_t integP[13][10];     // integral for phi (calculated from bin 1)
+  Double_t integPE[13][10];    // error
+  Double_t integPS[13][10];     // integral for phi (calculated from bin 1) (shuffled)
+  Double_t integPSE[13][10];    // error (shuffled)
+
+
+  while ( (key = (TKey*)nextkey())) {
+
+    list = (TList*)key->ReadObj();
+    listName = TString(list->GetName());
+
+    cout<<"Processing list "<<listName<<endl;
+
+    // ----------------------------------------------------
+    // plot QA histograms
+    if(listName.Contains("QA")){   
+
+      if(iList<13) iList++;
+      else{
+       cerr<<"TOO MANY LISTS!!!"<<endl;
+       return;
+      }  
+      
+      cQA[iList][iCanvas] = new TCanvas(listName,listName);
+      cQA[iList][iCanvas]->Divide(4,3);
+
+     cQA[iList][iCanvas]->cd(1);
+      TH1F* histEventStats = (TH1F*)list->FindObject("fHistEventStats");
+      if(histEventStats){
+       histEventStats->SetFillColor(9);
+       histEventStats->Draw();
+      }
+
+      cQA[iList][iCanvas]->cd(2);
+      TH1F* histTriggerStats = (TH1F*)list->FindObject("fHistTriggerStats");
+      if(histTriggerStats){
+       histTriggerStats->SetFillColor(9);
+       histTriggerStats->Draw();
+      }
+
+     cQA[iList][iCanvas]->cd(3);
+      TH1F* histTrackStats = (TH1F*)list->FindObject("fHistTrackStats");
+      if(histTrackStats){
+       histTrackStats->SetFillColor(9);
+       histTrackStats->Draw();
+      }
+
+      cQA[iList][iCanvas]->cd(4);
+      TH1F* histCentStats = (TH1F*)list->FindObject("fHistCentStats");
+      if(histCentStats){
+       histCentStats->SetFillColor(9);
+       histCentStats->Draw("colz");
+      }
+
+
+
+      cQA[iList][iCanvas]->cd(5);
+      TH1F* histVx = (TH1F*)list->FindObject("fHistVx");
+      if(histVx){
+       histVx->SetFillColor(9);
+       histVx->Draw();
+      }
+      cQA[iList][iCanvas]->cd(6);
+      TH1F* histVy = (TH1F*)list->FindObject("fHistVy");
+      if(histVy){
+       histVy->SetFillColor(9);
+       histVy->Draw();
+      }
+      
+      cQA[iList][iCanvas]->cd(7);
+      TH1F* histVz = (TH1F*)list->FindObject("fHistVz");
+      if(histVz){
+       histVz->SetFillColor(9);
+       histVz->Draw();
+      }
+
+      cQA[iList][iCanvas]->cd(8);
+      cQA[iList][iCanvas]->cd(8)->SetLogz();
+      TH2F* histDCA = (TH2F*)list->FindObject("fHistDCA");  
+      if(histDCA) histDCA->Draw("colz");
+      
+
+      cQA[iList][iCanvas]->cd(9);
+      cQA[iList][iCanvas]->cd(9)->SetLogz();
+      TH2F* histClus = (TH2F*)list->FindObject("fHistClus");
+      if(histClus) histClus->Draw("colz");
+      
+      cQA[iList][iCanvas]->cd(10);
+      TH1F* histPt = (TH1F*)list->FindObject("fHistPt");
+      if(histPt){
+       histPt->SetFillColor(9);
+       histPt->Draw();
+      }
+      
+      cQA[iList][iCanvas]->cd(11);
+      TH1F* histEta = (TH1F*)list->FindObject("fHistEta");
+      if(histEta){
+       histEta->SetFillColor(9);
+       histEta->Draw();
+      }
+      
+      cQA[iList][iCanvas]->cd(12);
+      TH1F* histPhi = (TH1F*)list->FindObject("fHistPhi");
+      if(histPhi){
+       histPhi->SetFillColor(9);
+       histPhi->Draw();
+      }
+
+      // centrality estimator QA
+      cQAV0M->cd(iCanvas+1);
+      cQAV0M->cd(iCanvas+1)->SetLogz();
+      TH1F* histV0M = (TH1F*)list->FindObject("fHistV0M");
+      if(histV0M){
+       histV0M->Draw("colz");
+      }
+
+      cQARef->cd(iCanvas+1);
+      cQARef->cd(iCanvas+1)->SetLogz();
+      TH1F* histRef = (TH1F*)list->FindObject("fHistRefTracks");
+      if(histRef){
+       histRef->Draw("colz");
+      }
+    }
+    // ----------------------------------------------------
+
+    // ----------------------------------------------------
+    // calculate and plot BF 
+    if(listName.Contains("BF_")){
+
+      for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
+       
+       cBF[iList][iCanvas] = new TCanvas(Form("cBF_%d_%d",iList,iCanvas),Form("cBF_%d_%d",iList,iCanvas));
+       cBF[iList][iCanvas]->Divide(3,3);
+       
+      }
+      
+
+      for(Int_t a = 0; a < 7; a++){
+
+       cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;
+
+       // create the BF object
+       bf[iList][a]  = new AliBalance();
+
+       fHistP[iList][a] = (TH2D*)list->FindObject(Form("fHistP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistN[iList][a] = (TH2D*)list->FindObject(Form("fHistN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistPP[iList][a] = (TH2D*)list->FindObject(Form("fHistPP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistPN[iList][a] = (TH2D*)list->FindObject(Form("fHistPN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistNP[iList][a] = (TH2D*)list->FindObject(Form("fHistNP%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistNN[iList][a] = (TH2D*)list->FindObject(Form("fHistNN%s%s",gBFAnalysisType[a].Data(),centEst.Data()));
+
+       // rebin histograms (be careful with divider!)
+       if(a==6){
+         fHistP[iList][a]->RebinY(5);
+         fHistN[iList][a]->RebinY(5);
+         fHistPP[iList][a]->RebinY(5);
+         fHistPN[iList][a]->RebinY(5);
+         fHistNP[iList][a]->RebinY(5);
+         fHistNN[iList][a]->RebinY(5);
+       }
+       else{
+         fHistP[iList][a]->RebinY(fRebin);
+         fHistN[iList][a]->RebinY(fRebin);
+         fHistPP[iList][a]->RebinY(fRebin);
+         fHistPN[iList][a]->RebinY(fRebin);
+         fHistNP[iList][a]->RebinY(fRebin);
+         fHistNN[iList][a]->RebinY(fRebin);
+       }
+
+       fHistP[iList][a]->SetName(Form("%s_%d_%d",fHistP[iList][a]->GetName(),iList,iCanvas));
+       fHistN[iList][a]->SetName(Form("%s_%d_%d",fHistN[iList][a]->GetName(),iList,iCanvas));
+       fHistPP[iList][a]->SetName(Form("%s_%d_%d",fHistPP[iList][a]->GetName(),iList,iCanvas));
+       fHistPN[iList][a]->SetName(Form("%s_%d_%d",fHistPN[iList][a]->GetName(),iList,iCanvas));
+       fHistNP[iList][a]->SetName(Form("%s_%d_%d",fHistNP[iList][a]->GetName(),iList,iCanvas));
+       fHistNN[iList][a]->SetName(Form("%s_%d_%d",fHistNN[iList][a]->GetName(),iList,iCanvas));
+
+       // set histograms in AliBalance object
+       bf[iList][a]->SetHistNp(a, fHistP[iList][a]);
+       bf[iList][a]->SetHistNn(a, fHistN[iList][a]);
+       bf[iList][a]->SetHistNpp(a, fHistPP[iList][a]);
+       bf[iList][a]->SetHistNpn(a, fHistPN[iList][a]);
+       bf[iList][a]->SetHistNnp(a, fHistNP[iList][a]);
+       bf[iList][a]->SetHistNnn(a, fHistNN[iList][a]);
+
+       for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
+
+         gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow);
+         gbf[iList][iCanvas][a]->SetName(Form("BF_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
+
+         cBF[iList][iCanvas]->cd(a+1);
+         gbf[iList][iCanvas][a]->SetMarkerStyle(20);
+
+         if(!bHistos){
+           gbf[iList][iCanvas][a]->DrawCopy("AP");
+           if(a==1){
+             GetWeightedMean(gbf[iList][iCanvas][a],fStartBinBFWidth,WM[iList][iCanvas],WME[iList][iCanvas]); // for eta recalculate width (from 0.1 only!)
+             GetIntegral(gbf[iList][iCanvas][a],integ[iList][iCanvas],integE[iList][iCanvas]);
+           }
+           else if(a==6){
+             GetWeightedMean(gbf[iList][iCanvas][a],1,WMP[iList][iCanvas],WMPE[iList][iCanvas]); // for phi calculate width 
+             GetIntegral(gbf[iList][iCanvas][a],integP[iList][iCanvas],integPE[iList][iCanvas]);
+           }
+         }
+         else{
+           fHistPN[iList][a]->SetLineColor(2);
+           fHistPN[iList][a]->ProjectionY(Form("pn%d",a))->DrawCopy();
+           fHistPP[iList][a]->SetLineColor(1);
+           fHistPP[iList][a]->ProjectionY(Form("pp%d",a))->DrawCopy("same");
+           fHistNP[iList][a]->SetLineColor(4);
+           fHistNP[iList][a]->ProjectionY(Form("np%d",a))->DrawCopy("same");
+           fHistNN[iList][a]->SetLineColor(8);
+           fHistNN[iList][a]->ProjectionY(Form("nn%d",a))->DrawCopy("same");
+         }
+       }
+      }
+    }
+    // ----------------------------------------------------
+
+    // ----------------------------------------------------
+    // calculate and plot BF (shuffled)
+    if(listName.Contains("BFShuffled")&&listName.Contains(extraString.Data())){
+
+      for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
+       
+       cBFS[iList][iCanvas] = new TCanvas(Form("Shuffled_%d_%d",iList,iCanvas),Form("Shuffled_%d_%d",iList,iCanvas));
+       cBFS[iList][iCanvas]->Divide(3,3);
+       
+      }
+
+      for(Int_t a = 0; a < 7; a++){
+
+       // create the BF object
+       bfs[iList][a]  = new AliBalance();
+
+       fHistPS[iList][a] = (TH2D*)list->FindObject(Form("fHistP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistNS[iList][a] = (TH2D*)list->FindObject(Form("fHistN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistPPS[iList][a] = (TH2D*)list->FindObject(Form("fHistPP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistPNS[iList][a] = (TH2D*)list->FindObject(Form("fHistPN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistNPS[iList][a] = (TH2D*)list->FindObject(Form("fHistNP%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
+       fHistNNS[iList][a] = (TH2D*)list->FindObject(Form("fHistNN%s_shuffle%s",gBFAnalysisType[a].Data(),centEst.Data()));
+
+       // rebin histograms (be careful with divider!)
+       if(a==6){
+         fHistPS[iList][a]->RebinY(5);
+         fHistNS[iList][a]->RebinY(5);
+         fHistPPS[iList][a]->RebinY(5);
+         fHistPNS[iList][a]->RebinY(5);
+         fHistNPS[iList][a]->RebinY(5);
+         fHistNNS[iList][a]->RebinY(5);
+       }
+       else{
+         fHistPS[iList][a]->RebinY(fRebin);
+         fHistNS[iList][a]->RebinY(fRebin);
+         fHistPPS[iList][a]->RebinY(fRebin);
+         fHistPNS[iList][a]->RebinY(fRebin);
+         fHistNPS[iList][a]->RebinY(fRebin);
+         fHistNNS[iList][a]->RebinY(fRebin);
+       }
+
+       fHistPS[iList][a]->SetName(Form("%s_%d_%d",fHistPS[iList][a]->GetName(),iList,iCanvas));
+       fHistNS[iList][a]->SetName(Form("%s_%d_%d",fHistNS[iList][a]->GetName(),iList,iCanvas));
+       fHistPPS[iList][a]->SetName(Form("%s_%d_%d",fHistPPS[iList][a]->GetName(),iList,iCanvas));
+       fHistPNS[iList][a]->SetName(Form("%s_%d_%d",fHistPNS[iList][a]->GetName(),iList,iCanvas));
+       fHistNPS[iList][a]->SetName(Form("%s_%d_%d",fHistNPS[iList][a]->GetName(),iList,iCanvas));
+       fHistNNS[iList][a]->SetName(Form("%s_%d_%d",fHistNNS[iList][a]->GetName(),iList,iCanvas));
+
+       // set histograms in AliBalance object
+       bfs[iList][a]->SetHistNp(a, fHistPS[iList][a]);
+       bfs[iList][a]->SetHistNn(a, fHistNS[iList][a]);
+       bfs[iList][a]->SetHistNpp(a, fHistPPS[iList][a]);
+       bfs[iList][a]->SetHistNpn(a, fHistPNS[iList][a]);
+       bfs[iList][a]->SetHistNnp(a, fHistNPS[iList][a]);
+       bfs[iList][a]->SetHistNnn(a, fHistNNS[iList][a]);
+
+       for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
+         
+         gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow);
+         gbfs[iList][iCanvas][a]->SetName(Form("BFS_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
+         
+         cBFS[iList][iCanvas]->cd(a+1);
+         gbfs[iList][iCanvas][a]->SetMarkerStyle(20);
+         if(!bHistos){
+           gbfs[iList][iCanvas][a]->DrawCopy("AP");
+           if(a==1){
+             GetWeightedMean(gbfs[iList][iCanvas][a],fStartBinBFWidth,WMS[iList][iCanvas],WMSE[iList][iCanvas]); 
+             GetIntegral(gbfs[iList][iCanvas][a],integS[iList][iCanvas],integSE[iList][iCanvas]); 
+           }
+           else if(a==6){
+             GetWeightedMean(gbfs[iList][iCanvas][a],1,WMPS[iList][iCanvas],WMPSE[iList][iCanvas]); // for phi calculate width 
+             GetIntegral(gbfs[iList][iCanvas][a],integPS[iList][iCanvas],integPSE[iList][iCanvas]); 
+           }
+         }
+         else{
+           fHistPNS[iList][a]->SetLineColor(2);
+           fHistPNS[iList][a]->ProjectionY(Form("pns%d",a))->DrawCopy();
+           fHistPPS[iList][a]->SetLineColor(1);
+           fHistPPS[iList][a]->ProjectionY(Form("pps%d",a))->DrawCopy("same");
+           fHistNPS[iList][a]->SetLineColor(4);
+           fHistNPS[iList][a]->ProjectionY(Form("nps%d",a))->DrawCopy("same");
+           fHistNNS[iList][a]->SetLineColor(8);
+           fHistNNS[iList][a]->ProjectionY(Form("nns%d",a))->DrawCopy("same");
+         }
+       }
+      }
+    }
+    // ----------------------------------------------------
+  }
+  
+  // for BF calculation create also graphs with weighted mean for eta
+  if(iCanvas == 0) return; 
+
+  TGraphErrors *gWM[13];
+  TGraphErrors *gWMS[13];
+  TGraphErrors *gWMP[13];
+  TGraphErrors *gWMPS[13];
+  TGraphErrors *ginteg[13];
+  TGraphErrors *gintegS[13];
+  TGraphErrors *gintegP[13];
+  TGraphErrors *gintegPS[13];
+  for(Int_t i = 0; i < 13; i++){
+    gWM[i] = NULL;
+    gWMS[i] = NULL;
+    ginteg[i] = NULL;
+    gintegS[i] = NULL;
+    gWMP[i] = NULL;
+    gWMPS[i] = NULL;
+    gintegP[i] = NULL;
+    gintegPS[i] = NULL;
+  }
+
+
+  if(!bHistos){
+    for(Int_t i = 0; i < iList+1; i++){
+      gWM[i] = new TGraphErrors(iCanvas,cent,WM[i],centE,WME[i]);
+      if(iList==1) gWM[i]->SetName("gCentrality");  
+      else gWM[i]->SetName(Form("gCentrality_%d",i));
+      gWMS[i] = new TGraphErrors(iCanvas,cent,WMS[i],centE,WMSE[i]); 
+      if(iList==1) gWMS[i]->SetName("gCentralityS");  
+      else gWMS[i]->SetName(Form("gCentralityS_%d",i)); 
+
+      gWMP[i] = new TGraphErrors(iCanvas,cent,WMP[i],centE,WMPE[i]);
+      if(iList==1) gWMP[i]->SetName("gCentralityPhi");  
+      else gWMP[i]->SetName(Form("gCentralityPhi_%d",i));
+      gWMPS[i] = new TGraphErrors(iCanvas,cent,WMPS[i],centE,WMPSE[i]); 
+      if(iList==1) gWMPS[i]->SetName("gCentralityPhiS");  
+      else gWMPS[i]->SetName(Form("gCentralityPhiS_%d",i)); 
+
+      ginteg[i] = new TGraphErrors(iCanvas,cent,integ[i],centE,integE[i]);
+      if(iList==1) ginteg[i]->SetName("gIntegral");  
+      else ginteg[i]->SetName(Form("gIntegral_%d",i));
+      gintegS[i] = new TGraphErrors(iCanvas,cent,integS[i],centE,integSE[i]); 
+      if(iList==1) gintegS[i]->SetName("gIntegralS");  
+      else gintegS[i]->SetName(Form("gIntegralS_%d",i)); 
+
+      gintegP[i] = new TGraphErrors(iCanvas,cent,integP[i],centE,integPE[i]);
+      if(iList==1) gintegP[i]->SetName("gIntegraPhil");  
+      else gintegP[i]->SetName(Form("gIntegralPhi_%d",i));
+      gintegPS[i] = new TGraphErrors(iCanvas,cent,integPS[i],centE,integPSE[i]); 
+      if(iList==1) gintegPS[i]->SetName("gIntegralPhiS");  
+      else gintegPS[i]->SetName(Form("gIntegralPhiS_%d",i)); 
+    }
+  }
+
+  TFile *fOut = TFile::Open(Form("Histograms_WMstart%d_rebin%d_%s", fStartBinBFWidth, fRebin,inFile.Data()),"RECREATE");
+  fOut->cd();
+  for(Int_t i = 0; i < iList+1; i++){
+    for(Int_t a = 0; a < 7; a++){
+      
+      if(fHistPN[i][a]){
+       (fHistPN[i][a]->ProjectionY(Form("hPN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistPP[i][a]->ProjectionY(Form("hPP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistNP[i][a]->ProjectionY(Form("hNP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistNN[i][a]->ProjectionY(Form("hNN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistP[i][a]->ProjectionY(Form("hP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistN[i][a]->ProjectionY(Form("hN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+      }
+
+      if(fHistPNS[i][a]){
+       (fHistPNS[i][a]->ProjectionY(Form("hPNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistPPS[i][a]->ProjectionY(Form("hPPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistNPS[i][a]->ProjectionY(Form("hNPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistNNS[i][a]->ProjectionY(Form("hNNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistPS[i][a]->ProjectionY(Form("hPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+       (fHistNS[i][a]->ProjectionY(Form("hNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
+      }
+      
+      for(Int_t j = 0; j < iCanvas; j++){
+       
+       if(gbf[i][j][a]){
+         //printout in text format for delta eta
+         if(a==1){
+           cout<<"Balance Function "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
+           cout<<endl;
+           cout<<"Balance Function error "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
+           cout<<endl;
+         } 
+         else if(a==6){
+           cout<<"Balance Function Phi "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
+           cout<<endl;
+           cout<<"Balance Function Phi error "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX();k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
+           cout<<endl;
+         } 
+         gbf[i][j][a]->Write();
+         gbf[i][j][a]->Delete();
+       }
+       if(gbfs[i][j][a]){
+         if(a==1){
+           cout<<"Balance Function (shuffled) "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
+           cout<<endl;
+           cout<<"Balance Function error (shuffled) "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
+           cout<<endl;
+         } 
+         else if(a==6){
+           cout<<"Balance Function Phi (shuffled) "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
+           cout<<endl;
+           cout<<"Balance Function Phi error (shuffled) "<<i<<" "<<j<<" : ";
+           for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX();k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
+           cout<<endl;
+         } 
+         gbfs[i][j][a]->Write();
+         gbfs[i][j][a]->Delete();
+       }
+      }
+    }
+
+    Double_t x,y;
+      
+    if(gWM[i]){
+      cout<<"Balance Function WM "<<i<<" : ";
+      for(Int_t k = 0; k < gWM[i]->GetN();k++){
+       gWM[i]->GetPoint(k,x,y);
+       cout<<y<<", ";      
+      }
+      cout<<endl;
+      cout<<"Balance Function WM error "<<i<<" : ";
+      for(Int_t k = 0; k < gWM[i]->GetN();k++){
+       cout<<gWM[i]->GetErrorY(k)<<", ";      
+      }
+      cout<<endl;
+      gWM[i]->Write();
+    } 
+    if(gWMS[i]){
+      cout<<"Balance Function WM (shuffled) "<<i<<" : ";
+      for(Int_t k = 0; k < gWMS[i]->GetN();k++){
+       gWMS[i]->GetPoint(k,x,y);
+       cout<<y<<", ";      
+      }
+      cout<<endl;
+      cout<<"Balance Function WM error (shuffled) "<<i<<" : ";
+      for(Int_t k = 0; k < gWMS[i]->GetN();k++){
+       cout<<gWMS[i]->GetErrorY(k)<<", ";      
+      }
+      cout<<endl;
+      gWMS[i]->Write();
+    } 
+
+   if(gWMP[i]){
+      cout<<"Balance Function Phi WMP "<<i<<" : ";
+      for(Int_t k = 0; k < gWMP[i]->GetN();k++){
+       gWMP[i]->GetPoint(k,x,y);
+       cout<<y<<", ";      
+      }
+      cout<<endl;
+      cout<<"Balance Function Phi WMP error "<<i<<" : ";
+      for(Int_t k = 0; k < gWMP[i]->GetN();k++){
+       cout<<gWMP[i]->GetErrorY(k)<<", ";      
+      }
+      cout<<endl;
+      gWMP[i]->Write();
+    } 
+    if(gWMPS[i]){
+      cout<<"Balance Function Phi WMP (shuffled) "<<i<<" : ";
+      for(Int_t k = 0; k < gWMPS[i]->GetN();k++){
+       gWMPS[i]->GetPoint(k,x,y);
+       cout<<y<<", ";      
+      }
+      cout<<endl;
+      cout<<"Balance Function Phi WMP error (shuffled) "<<i<<" : ";
+      for(Int_t k = 0; k < gWMPS[i]->GetN();k++){
+       cout<<gWMPS[i]->GetErrorY(k)<<", ";      
+      }
+      cout<<endl;
+      gWMPS[i]->Write();
+    } 
+
+    if(ginteg[i]) ginteg[i]->Write();
+    if(gintegS[i]) gintegS[i]->Write();
+    if(gintegP[i]) gintegP[i]->Write();
+    if(gintegPS[i]) gintegPS[i]->Write();
+  }
+  fOut->Close();
+  f->Close();
+  gROOT->Reset();
+}
+
+//____________________________________________________________________//
+void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Double_t &WM, Double_t &WME) {
+
+  //Prints the calculated width of the BF and its error
+  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
+  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
+  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
+  Double_t deltaBalP2 = 0.0, integral = 0.0;
+  Double_t deltaErrorNew = 0.0;
+
+  //Retrieve this variables from Histogram
+  Int_t fNumberOfBins = gHistBalance->GetNbinsX();
+  Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
+  
+  cout<<"=================================================="<<endl;
+  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
+  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
+  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
+    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
+  } 
+  cout<<"=================================================="<<endl;
+  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
+    gSumXi += gHistBalance->GetBinCenter(i);
+    gSumBi += gHistBalance->GetBinContent(i);
+    gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
+    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
+    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
+    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
+    gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
+    
+    deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
+    integral += fP2Step*gHistBalance->GetBinContent(i);
+  }
+  for(Int_t i = fStartBin; i < fNumberOfBins; i++)
+    deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
+  
+  Double_t integralError = TMath::Sqrt(deltaBalP2);
+  
+  Double_t delta = gSumBiXi / gSumBi;
+  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
+  
+  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+  cout<<"New error: "<<deltaErrorNew<<endl;
+  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+  cout<<"=================================================="<<endl;
+
+  WM  = delta;
+  WME = deltaError;
+}
+
+//____________________________________________________________________//
+void GetIntegral(TH1D *gHistBalance, Double_t &integ, Double_t &integE) {
+
+  //always start with 1
+  Int_t fStartBin = 1; 
+
+  //Prints the calculated width of the BF and its error
+  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
+  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
+  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
+  Double_t deltaBalP2 = 0.0, integral = 0.0;
+  Double_t deltaErrorNew = 0.0;
+
+  //Retrieve this variables from Histogram
+  Int_t fNumberOfBins = gHistBalance->GetNbinsX();
+  Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
+  
+  cout<<"=================================================="<<endl;
+  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
+  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
+  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
+    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
+  } 
+  cout<<"=================================================="<<endl;
+  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
+    gSumXi += gHistBalance->GetBinCenter(i);
+    gSumBi += gHistBalance->GetBinContent(i);
+    gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
+    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
+    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
+    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
+    gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
+    
+    deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
+    integral += fP2Step*gHistBalance->GetBinContent(i);
+  }
+  for(Int_t i = fStartBin; i < fNumberOfBins; i++)
+    deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
+  
+  Double_t integralError = TMath::Sqrt(deltaBalP2);
+  
+  Double_t delta = gSumBiXi / gSumBi;
+  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
+  
+  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+  cout<<"New error: "<<deltaErrorNew<<endl;
+  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+  cout<<"=================================================="<<endl;
+
+  integ  = integral;
+  integE = integralError;
+}
+
+//___________________________________________________________//
+// NOT USED any more
+//___________________________________________________________//
+void mergeOutput(const char* outputDir) {
+  //Function to merge the output of the sub-jobs
+  //Create a BF object
+  AliBalance *bf = new AliBalance();
+
+  //connect to AliEn's API services
+  TGrid::Connect("alien://"); 
+
+  //Getting the output dir from the env. variable 
+  //(JDL field: JDLVariables={"OutputDir"};)
+  TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");
+  
+  Int_t nEntries = result->GetEntries();
+
+  TString alienUrl;
+  TDirectoryFile *dirSubJob;
+
+  TString gCutName[4] = {"Total","Offline trigger",
+                         "Vertex","Analyzed"};
+  TH1F *fHistEventStats = new TH1F("fHistEventStats",
+                                  "Event statistics;;N_{events}",
+                                  4,0.5,4.5);
+  for(Int_t i = 1; i <= 4; i++)
+    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
+
+  AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");
+  for(Int_t i = 0; i < nEntries; i++) {
+    alienUrl = result->GetKey(i,"turl");
+    alienUrl += "#AnalysisResults.root";
+    Printf("Opening file: %s",alienUrl.Data());
+    TFile *file = TFile::Open(alienUrl.Data());
+    dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWGCFEbyE.outputBalanceFunctionAnalysis.root"));
+
+    //merge BF
+    AliBalance *bfSubJob = dynamic_cast<AliBalance *>(dirSubJob->Get("AliBalance"));
+    //bfSubJob->PrintResults();
+    bf->Merge(bfSubJob);
+    //delete bfSubJob;
+
+    //merge event stats
+    TList *listSubJob = dynamic_cast<TList *>(dirSubJob->Get("listQA"));
+    fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));
+
+    bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));
+    delete listSubJob;
+  }
+
+  //Create the output file
+  TString outputFile = "AnalysisResults.Merged.root";
+  TFile *foutput = TFile::Open(outputFile.Data(),"recreate");
+  TDirectoryFile *dirOutput = new TDirectoryFile();
+  dirOutput->SetName("PWGCFEbyE.outputBalanceFunctionAnalysis.root");
+  //dirOutput->cd();
+  dirOutput->Add(bf);
+  TList *list = new TList();
+  list->SetName("listQA");
+  list->Add(fHistEventStats);
+  list->Add(bfTrackCuts);
+  dirOutput->Add(list);
+  dirOutput->Write();
+  bf->Write();
+  list->Write();
+  foutput->Close();
+
+    //cout<<alienUrl.Data()<<endl;
+}