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