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