PWGCF/EBYE/macros/readBalanceFunction.C

   1 const TString gBFAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
   2
   3 const Int_t nrOfCentralities = 9;
   4 const Double_t centralityArray[nrOfCentralities+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};  // in centrality percentile (0-5,5-10,10-20,20-30,...,70-80)
   5 //const Double_t centralityArray[nrOfCentralities+1] = {0.,1.,2.,3.,4.,6.,10.,20.,30.,80.};  // in centrality percentile (0-5,5-10,10-20,20-30,...,70-80)
   6 const Double_t cent[nrOfCentralities]  = {382.8,329.7,260.5,186.4,128.9,85.,52.8,30.,15.8};   // hard coded at the moment for centrality percentiles
   7 const Double_t centE[nrOfCentralities] = {3.1,4.6,4.4,3.9,3.3,2.6,2.0,1.3,0.6};               // (0-5,5-10,10-20,20-30,...,70-80)
   8
   9 void readBalanceFunction(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root",Int_t fStartBinBFWidth = 3, Int_t fRebin = 2,Int_t fStartBinBFWidthPhi = 2, Int_t fRebinPhi = 2,TString centEst = "V0M",Double_t etaWindow = -1, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE) {
  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,etaBins, correctWithEfficiency,correctWithAcceptanceOnly);
  24
  25
  26   //Merge the output
  27   //mergeOutput("/alice/cern.ch/user/p/pchrist/Balance/pp/7TeV/LHC10b/output/");
  28 }
  29
  30 //___________________________________________________________//
  31 void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", Int_t fStartBinBFWidth = 1, Int_t fRebin = 1, Int_t fStartBinBFWidthPhi = 1, Int_t fRebinPhi = 1, TString centEst = "V0M",TString extraString = "", Double_t etaWindow = -1, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE) {
  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,correctWithEfficiency,correctWithAcceptanceOnly);
 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,etaBins,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,-1,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, correctWithEfficiency,correctWithAcceptanceOnly);
 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,etaBins,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,-1,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 = NULL;
 535   if(etaWindow > 0 || etaBins > -1){
 536     if(correctWithEfficiency){
 537       if(correctWithAcceptanceOnly){
 538         fOut = TFile::Open(Form("Histograms_AccCorrWithAcceptanceOnly_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
 539       }
 540       else{
 541         fOut = TFile::Open(Form("Histograms_AccCorrWithEfficiency_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
 542       }
 543     }
 544     else{
 545       fOut = TFile::Open(Form("Histograms_AccCorr_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
 546     }
 547   }
 548   else{
 549     fOut = TFile::Open(Form("Histograms_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
 550   }
 551   fOut->cd();
 552   for(Int_t i = 0; i < iList+1; i++){
 553     cout<<"PROCESS LIST "<<i<<" NOW!"<<endl;
 554     for(Int_t a = 0; a < 7; a++){
 555     cout<<"PROCESS VARIABLE "<<a<<" NOW!"<<endl;
 556
 557       if(fHistPN[i][a]){
 558         (fHistPN[i][a]->ProjectionY(Form("hPN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 559         (fHistPP[i][a]->ProjectionY(Form("hPP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 560         (fHistNP[i][a]->ProjectionY(Form("hNP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 561         (fHistNN[i][a]->ProjectionY(Form("hNN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 562         (fHistP[i][a]->ProjectionY(Form("hP_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 563         (fHistN[i][a]->ProjectionY(Form("hN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 564       }
 565
 566       if(fHistPNS[i][a]){
 567         (fHistPNS[i][a]->ProjectionY(Form("hPNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 568         (fHistPPS[i][a]->ProjectionY(Form("hPPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 569         (fHistNPS[i][a]->ProjectionY(Form("hNPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 570         (fHistNNS[i][a]->ProjectionY(Form("hNNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 571         (fHistPS[i][a]->ProjectionY(Form("hPS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 572         (fHistNS[i][a]->ProjectionY(Form("hNS_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
 573       }
 574
 575       //printout in text format for delta eta
 576       for(Int_t j = 0; j < iCanvas; j++){
 577         cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
 578
 579         if(gbf[i][j][a]){
 580           if(a==1){
 581             cout<<"Double_t gALICEDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaEta] = {";
 582             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
 583             cout<<gbf[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
 584             cout<<"Double_t gALICEDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaEta] = {";
 585             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
 586             cout<<gbf[i][j][a]->GetBinError(k+1)<<"};"<<endl;
 587           }
 588           else if(a==6){
 589             cout<<"Double_t gALICEDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaPhi] = {";
 590             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinContent(k+1)<<", ";
 591             cout<<gbf[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
 592             cout<<"Double_t gALICEDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaPhi] = {";
 593             for(Int_t k = 0; k < gbf[i][j][a]->GetNbinsX()-1;k++) cout<<gbf[i][j][a]->GetBinError(k+1)<<", ";
 594             cout<<gbf[i][j][a]->GetBinError(k+1)<<"};"<<endl;
 595           }
 596           gbf[i][j][a]->Write();
 597           gbf[i][j][a]->Delete();
 598         }
 599         if(gbfs[i][j][a]){
 600           if(a==1){
 601             cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaEta] = {";
 602             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
 603             cout<<gbfs[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
 604             cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaEtaCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaEta] = {";
 605             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
 606             cout<<gbfs[i][j][a]->GetBinError(k+1)<<"};"<<endl;
 607           }
 608           else if(a==6){
 609             cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"[nBinsInDeltaPhi] = {";
 610             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinContent(k+1)<<", ";
 611             cout<<gbfs[i][j][a]->GetBinContent(k+1)<<"};"<<endl;
 612             cout<<"Double_t gALICEShuffledDataBalanceFunctionInDeltaPhiCentrality"<<centralityArray[j]<<"to"<<centralityArray[j+1]<<"Error[nBinsInDeltaPhi] = {";
 613             for(Int_t k = 0; k < gbfs[i][j][a]->GetNbinsX()-1;k++) cout<<gbfs[i][j][a]->GetBinError(k+1)<<", ";
 614             cout<<gbfs[i][j][a]->GetBinError(k+1)<<"};"<<endl;
 615           }
 616           gbfs[i][j][a]->Write();
 617           gbfs[i][j][a]->Delete();
 618         }
 619         cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
 620         cout<<endl;
 621       }
 622     }
 623
 624     Double_t x,y;
 625
 626     cout<<"//================================ALICE================================//"<<endl;
 627     if(gWM[i]){
 628       cout<<"Double_t gWeightedMeanInEtaAlice[nCentralityBins] = {";
 629       for(Int_t k = 0; k < gWM[i]->GetN()-1;k++){
 630         gWM[i]->GetPoint(k,x,y);
 631         cout<<y<<", ";
 632       }
 633       gWM[i]->GetPoint(k,x,y);
 634       cout<<y<<"};"<<endl;
 635
 636       cout<<"Double_t gWeightedMeanInEtaAliceError[nCentralityBins] = {";
 637       for(Int_t k = 0; k < gWM[i]->GetN()-1;k++){
 638         cout<<gWM[i]->GetErrorY(k)<<", ";
 639       }
 640       cout<<gWM[i]->GetErrorY(k)<<"};"<<endl;
 641       gWM[i]->Write();
 642     }
 643     if(gWMS[i]){
 644       cout<<"Double_t gShuffledWeightedMeanInEtaAlice[nCentralityBins] = {";
 645       for(Int_t k = 0; k < gWMS[i]->GetN()-1;k++){
 646         gWMS[i]->GetPoint(k,x,y);
 647         cout<<y<<", ";
 648       }
 649       gWMS[i]->GetPoint(k,x,y);
 650       cout<<y<<"};"<<endl;
 651
 652       cout<<"Double_t gShuffledWeightedMeanInEtaAliceError[nCentralityBins] = {";
 653       for(Int_t k = 0; k < gWMS[i]->GetN()-1;k++){
 654         cout<<gWMS[i]->GetErrorY(k)<<", ";
 655       }
 656       cout<<gWMS[i]->GetErrorY(k)<<"};"<<endl;
 657       cout<<endl;
 658       gWMS[i]->Write();
 659     }
 660
 661     if(gWMP[i]){
 662       cout<<"Double_t gWeightedMeanInPhiAlice[nCentralityBins] = {";
 663       for(Int_t k = 0; k < gWMP[i]->GetN()-1;k++){
 664         gWMP[i]->GetPoint(k,x,y);
 665         cout<<y<<", ";
 666       }
 667       gWMP[i]->GetPoint(k,x,y);
 668       cout<<y<<"};"<<endl;
 669
 670       cout<<"Double_t gWeightedMeanInPhiAliceError[nCentralityBins] = {";
 671       for(Int_t k = 0; k < gWMP[i]->GetN()-1;k++){
 672         cout<<gWMP[i]->GetErrorY(k)<<", ";
 673       }
 674       cout<<gWMP[i]->GetErrorY(k)<<"};"<<endl;
 675       gWMP[i]->Write();
 676     }
 677     if(gWMPS[i]){
 678       cout<<"Double_t gShuffledWeightedMeanInPhiAlice[nCentralityBins] = {";
 679       for(Int_t k = 0; k < gWMPS[i]->GetN()-1;k++){
 680         gWMPS[i]->GetPoint(k,x,y);
 681         cout<<y<<", ";
 682       }
 683       gWMPS[i]->GetPoint(k,x,y);
 684       cout<<y<<"};"<<endl;
 685
 686       cout<<"Double_t gShuffledWeightedMeanInPhiAliceError[nCentralityBins] = {";
 687       for(Int_t k = 0; k < gWMPS[i]->GetN()-1;k++){
 688         cout<<gWMPS[i]->GetErrorY(k)<<", ";
 689       }
 690       cout<<gWMPS[i]->GetErrorY(k)<<"};"<<endl;
 691       cout<<endl;
 692       gWMPS[i]->Write();
 693     }
 694     if(ginteg[i]){
 695       cout<<"Double_t gIntegralEtaAlice[nCentralityBins] = {";
 696       for(Int_t k = 0; k < ginteg[i]->GetN()-1;k++){
 697         ginteg[i]->GetPoint(k,x,y);
 698         cout<<y<<", ";
 699       }
 700       ginteg[i]->GetPoint(k,x,y);
 701       cout<<y<<"};"<<endl;
 702       cout<<endl;
 703
 704       cout<<"Double_t gIntegralErrorEtaAlice[nCentralityBins] = {";
 705       for(Int_t k = 0; k < ginteg[i]->GetN()-1;k++){
 706         cout<<ginteg[i]->GetErrorY(k)<<", ";
 707       }
 708       cout<<ginteg[i]->GetErrorY(k)<<"};"<<endl;
 709       cout<<endl;
 710
 711       ginteg[i]->Write();
 712
 713       cout<<"Double_t gIntegralPhiAlice[nCentralityBins] = {";
 714       for(Int_t k = 0; k < gintegP[i]->GetN()-1;k++){
 715         gintegP[i]->GetPoint(k,x,y);
 716         cout<<y<<", ";
 717       }
 718       gintegP[i]->GetPoint(k,x,y);
 719       cout<<y<<"};"<<endl;
 720       cout<<endl;
 721
 722       cout<<"Double_t gIntegralErrorPhiAlice[nCentralityBins] = {";
 723       for(Int_t k = 0; k < gintegP[i]->GetN()-1;k++){
 724         cout<<gintegP[i]->GetErrorY(k)<<", ";
 725       }
 726       cout<<gintegP[i]->GetErrorY(k)<<"};"<<endl;
 727       cout<<endl;
 728
 729       gintegP[i]->Write();
 730
 731       cout<<"Double_t gIntegralEtaShuffledAlice[nCentralityBins] = {";
 732       for(Int_t k = 0; k < gintegS[i]->GetN()-1;k++){
 733         gintegS[i]->GetPoint(k,x,y);
 734         cout<<y<<", ";
 735       }
 736       gintegS[i]->GetPoint(k,x,y);
 737       cout<<y<<"};"<<endl;
 738       cout<<endl;
 739
 740       cout<<"Double_t gIntegralErrorEtaShuffledAlice[nCentralityBins] = {";
 741       for(Int_t k = 0; k < gintegS[i]->GetN()-1;k++){
 742         cout<<gintegS[i]->GetErrorY(k)<<", ";
 743       }
 744       cout<<gintegS[i]->GetErrorY(k)<<"};"<<endl;
 745       cout<<endl;
 746
 747       gintegS[i]->Write();
 748
 749       cout<<"Double_t gIntegralPhiShuffledAlice[nCentralityBins] = {";
 750       for(Int_t k = 0; k < gintegPS[i]->GetN()-1;k++){
 751         gintegPS[i]->GetPoint(k,x,y);
 752         cout<<y<<", ";
 753       }
 754       gintegPS[i]->GetPoint(k,x,y);
 755       cout<<y<<"};"<<endl;
 756       cout<<endl;
 757
 758       cout<<"Double_t gIntegralErrorPhiShuffledAlice[nCentralityBins] = {";
 759       for(Int_t k = 0; k < gintegPS[i]->GetN()-1;k++){
 760         cout<<gintegPS[i]->GetErrorY(k)<<", ";
 761       }
 762       cout<<gintegPS[i]->GetErrorY(k)<<"};"<<endl;
 763       cout<<endl;
 764
 765       gintegPS[i]->Write();
 766
 767     }
 768   }
 769   fOut->Close();
 770   f->Close();
 771   gROOT->Reset();
 772 }
 773
 774 //____________________________________________________________________//
 775 void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME) {
 776
 777   //Prints the calculated width of the BF and its error
 778   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
 779   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
 780   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
 781   Double_t deltaBalP2 = 0.0, integral = 0.0;
 782   Double_t deltaErrorNew = 0.0;
 783
 784   //Retrieve this variables from Histogram
 785   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
 786   if(fStopBin > -1) fNumberOfBins = fStopBin;
 787   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
 788
 789   cout<<"=================================================="<<endl;
 790   cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
 791   cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
 792   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
 793     cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
 794   }
 795   cout<<"=================================================="<<endl;
 796   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
 797     gSumXi += gHistBalance->GetBinCenter(i);
 798     gSumBi += gHistBalance->GetBinContent(i);
 799     gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
 800     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 801     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 802     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
 803     gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
 804
 805     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
 806     integral += fP2Step*gHistBalance->GetBinContent(i);
 807   }
 808   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
 809     deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
 810
 811   Double_t integralError = TMath::Sqrt(deltaBalP2);
 812
 813   Double_t delta = gSumBiXi / gSumBi;
 814   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
 815
 816   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 817   cout<<"New error: "<<deltaErrorNew<<endl;
 818   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
 819   cout<<"=================================================="<<endl;
 820
 821   WM  = delta;
 822   WME = deltaError;
 823 }
 824
 825 //____________________________________________________________________//
 826 void GetIntegral(TH1D *gHistBalance, Double_t &integ, Double_t &integE) {
 827
 828   //always start with 1
 829   Int_t fStartBin = 1;
 830
 831   //Prints the calculated width of the BF and its error
 832   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
 833   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
 834   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
 835   Double_t deltaBalP2 = 0.0, integral = 0.0;
 836   Double_t deltaErrorNew = 0.0;
 837
 838   //Retrieve this variables from Histogram
 839   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
 840   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
 841
 842   cout<<"=================================================="<<endl;
 843   cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
 844   cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
 845   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
 846     cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
 847   }
 848   cout<<"=================================================="<<endl;
 849   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
 850     gSumXi += gHistBalance->GetBinCenter(i);
 851     gSumBi += gHistBalance->GetBinContent(i);
 852     gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
 853     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 854     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
 855     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
 856     gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
 857
 858     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
 859     integral += fP2Step*gHistBalance->GetBinContent(i);
 860   }
 861   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
 862     deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
 863
 864   Double_t integralError = TMath::Sqrt(deltaBalP2);
 865
 866   Double_t delta = gSumBiXi / gSumBi;
 867   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
 868
 869   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 870   cout<<"New error: "<<deltaErrorNew<<endl;
 871   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
 872   cout<<"=================================================="<<endl;
 873
 874   integ  = integral;
 875   integE = integralError;
 876 }
 877
 878 //___________________________________________________________//
 879 // NOT USED any more
 880 //___________________________________________________________//
 881 void mergeOutput(const char* outputDir) {
 882   //Function to merge the output of the sub-jobs
 883   //Create a BF object
 884   AliBalance *bf = new AliBalance();
 885
 886   //connect to AliEn's API services
 887   TGrid::Connect("alien://");
 888
 889   //Getting the output dir from the env. variable
 890   //(JDL field: JDLVariables={"OutputDir"};)
 891   TGridResult* result = gGrid->Query(outputDir,"*/root_archive.zip","","-l 1000");
 892
 893   Int_t nEntries = result->GetEntries();
 894
 895   TString alienUrl;
 896   TDirectoryFile *dirSubJob;
 897
 898   TString gCutName[4] = {"Total","Offline trigger",
 899                          "Vertex","Analyzed"};
 900   TH1F *fHistEventStats = new TH1F("fHistEventStats",
 901                                    "Event statistics;;N_{events}",
 902                                    4,0.5,4.5);
 903   for(Int_t i = 1; i <= 4; i++)
 904     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
 905
 906   AliESDtrackCuts *bfTrackCuts = new AliESDtrackCuts("bfTrackCuts");
 907   for(Int_t i = 0; i < nEntries; i++) {
 908     alienUrl = result->GetKey(i,"turl");
 909     alienUrl += "#AnalysisResults.root";
 910     Printf("Opening file: %s",alienUrl.Data());
 911     TFile *file = TFile::Open(alienUrl.Data());
 912     dirSubJob = dynamic_cast<TDirectoryFile *>(file->Get("PWGCFEbyE.outputBalanceFunctionAnalysis.root"));
 913
 914     //merge BF
 915     AliBalance *bfSubJob = dynamic_cast<AliBalance *>(dirSubJob->Get("AliBalance"));
 916     //bfSubJob->PrintResults();
 917     bf->Merge(bfSubJob);
 918     //delete bfSubJob;
 919
 920     //merge event stats
 921     TList *listSubJob = dynamic_cast<TList *>(dirSubJob->Get("listQA"));
 922     fHistEventStats->Add(dynamic_cast<TH1F *>(listSubJob->At(0)));
 923
 924     bfTrackCuts = dynamic_cast<AliESDtrackCuts *>(listSubJob->At(1));
 925     delete listSubJob;
 926   }
 927
 928   //Create the output file
 929   TString outputFile = "AnalysisResults.Merged.root";
 930   TFile *foutput = TFile::Open(outputFile.Data(),"recreate");
 931   TDirectoryFile *dirOutput = new TDirectoryFile();
 932   dirOutput->SetName("PWGCFEbyE.outputBalanceFunctionAnalysis.root");
 933   //dirOutput->cd();
 934   dirOutput->Add(bf);
 935   TList *list = new TList();
 936   list->SetName("listQA");
 937   list->Add(fHistEventStats);
 938   list->Add(bfTrackCuts);
 939   dirOutput->Add(list);
 940   dirOutput->Write();
 941   bf->Write();
 942   list->Write();
 943   foutput->Close();
 944
 945     //cout<<alienUrl.Data()<<endl;
 946 }