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