PWGCF/FLOW/macros/plotCumulants.C

   1 // Macro plotCumulants.C is used to represent true correlations, a.k.a. cumulants,
   2 // as a function of reference multiplicity. Different order cumulants provide
   3 // independent estimates of the flow harmonics. In the ideal case when only flow
   4 // correlations are present we have:
   5 //         c{2} = v^2, c{4} = -v^4, c{6} = 4v^6 and c{8} = -33v^8.
   6 // In practice cumulants can are calculated via Q-vectors (QC) or via formalism
   7 // of generating functions (GFC).
   8
   9 // Set how many output analysis files in total you want to access:
  10 const Int_t nFiles = 2;
  11
  12 // Set how many of those output analysis files you want to represent with a mesh (usually used to represent results of simulations):
  13 const Int_t nSim = 0;
  14
  15 // Set paths of all output analysis files (first the ones to be represented with mesh (simulations), then the ones to be represented with markers (real data))
  16 TString files[nFiles] = {"trackletsCorrectedNUA","default"}; // subdirectory names holding <commonOutputFileName>.root
  17 //TString files[nFiles] = {"outputCentrality0","outputCentrality1"}; // file names (use case: centrality train)
  18
  19 // Set analysis types for all output analysis files (can be "ESD","AOD","MC","","MK", ....):
  20 TString type[nFiles] = {"ESD","ESD"};
  21 //TString type[nFiles] = {"MK","MK"};
  22
  23 // Set mesh color:
  24 Int_t meshColor[nSim] = {};
  25
  26 // Set marker styles:
  27 Int_t markerStyle[nFiles-nSim] = {kFullSquare,kOpenSquare};
  28
  29 // Set marker colors:
  30 Int_t markerColor[nFiles-nSim] = {kBlack,kRed};
  31
  32 // Set legend entries:
  33 TString legendEntry[nFiles] = {"",""};
  34
  35 // Set if you want to rebin the histograms into wider multiplicity bins (set for each cumulant order separately):
  36 Bool_t rebin = kFALSE;
  37 Int_t nMergedBins[4] = {10,10,10,10}; // set how many original multiplicity bins will be merged into 1 new one
  38
  39 // Set if you whish to plot cumulants versus <reference multiplicity> (by default they are plotted versus # of RPs):
  40 Bool_t plotCumulantsVsReferenceMultiplicity = kFALSE;
  41 Bool_t showReferenceMultiplicityVsNoOfRPs = kFALSE;
  42
  43 // Set flow values whose theoretical contribution to cumulants will be shown on the plots with the straight coloured lines:
  44 Bool_t showTheoreticalLines = kFALSE;
  45 const Int_t nFlowValues = 1;
  46 Double_t v[nFlowValues] = {0.05};
  47 Int_t lineColor[nFlowValues] = {kRed};
  48
  49 // If the statistical error of 6th and 8th order cumulant is huge you may prefer not to show them:
  50 Bool_t plotOnly2ndAnd4thOrderCumulant = kFALSE;
  51
  52 // Set if you want independent canvas with results for reference flow vs multiplicity:
  53 Bool_t showRefFlowVsM = kTRUE;
  54 Int_t refFlowVsMMarkerStyle[nFiles] = {kFullSquare,kOpenSquare}; // marker style is different for different file
  55 Int_t refFlowVsMMarkerColor[4] = {kBlack,kRed,kBlue,kGreen+2}; // marker color is different for different cumulant order
  56 Int_t refFlowVsMeshOrder = -1; // set results of which order will be plotted as mesh in all pads [0=2nd, 1=4th, 2=6th, 3=8th, -1=do not show]
  57 Int_t refFlowVsMMeshColor = kRed-10; // mesh color for above specified order
  58 Double_t refFlowVsMxRange[2] = {1.,11000.}; // x range on the plots for reference multiplicity vs M
  59 Double_t refFlowVsMyRange[2] = {0.0,0.194}; // x range on the plots for reference multiplicity vs M
  60
  61 // Set if you want to show theoretical curves for the toy model:
  62 Bool_t showToyModel = kFALSE;
  63 const Int_t nToyModels = 2; // number of toy models with different parameters k, vn, v2n
  64 Double_t k[nToyModels] = {2.,4.};
  65 Double_t vn[nToyModels] = {0.25,0.25};
  66 Double_t v2n[nToyModels] = {0.0,0.0};
  67 Int_t lineColorToyModel[nToyModels] = {kBlack,kRed};
  68
  69 // For comparison sake show also GFC results with dotted line:
  70 Bool_t showAlsoGFCResults = kFALSE;
  71 Int_t gfcLineStyle = 3;
  72
  73 // For comparison sake show also Monte Carlo QC results with coloured mesh:
  74 Bool_t showAlsoMCEPResults = kFALSE;
  75 Bool_t showOnlyMCEPResults = kFALSE;
  76 Int_t mcepMeshColor[nSim] = {};
  77
  78 // Set the naming convention:
  79 Bool_t bLookInSubdirectories = kFALSE; // if kTRUE: Look in subdirectories <files[0]>, <files[1]>, ..., for output files <commonOutputFileName>
  80                                       // if kFALSE: Look in <pwd> for files <files[0]>.root, <files[1]>.root, ....
  81 TString commonOutputFileName = "AnalysisResults.root";
  82
  83 // Set method names which calculate cumulants vs multiplicity (do not touch these settings unless you are looking for a trouble):
  84 const Int_t nMethods = 3;
  85 TString method[nMethods] = {"QC","GFC","MCEP"};
  86
  87 TFile *commonOutputFiles[nFiles] = {NULL}; // common output files "AnalysisResults.root"
  88 TList *lists[nFiles][nMethods] = {{NULL}}; // lists cobj<method> holding objects with results for each method
  89 TH1D *cumulantsVsM[nFiles][nMethods][4] = {{{NULL}}}; // histograms with results for cumulants vs multiplicity (4 stands for 4 cumulant orders)
  90 TH1D *refFlowVsM[nFiles][nMethods][4] = {{{NULL}}}; // histograms with results for reference flow vs multiplicity (4 stands for 4 cumulant orders)
  91 TGraph *lines[nFlowValues][4] = {{NULL}}; // lines denoting theoretical flow contribution to cumulants
  92 TProfile *refMultVsNoOfRPs[nFiles] = {NULL}; // <reference multipicity> versus # of RPs
  93 TF1 *toyModel[2][nToyModels] = {{NULL}}; // [cumulant order][number of toy models]
  94
  95 // Ranges for plots:
  96 Double_t xMin[4]={0.};
  97 Double_t xMax[4]={0.};
  98 Double_t yMin[4]={0.};
  99 Double_t yMax[4]={0.};
 100
 101 enum libModes {mLocal,mLocalSource};
 102
 103 void plotCumulants(Int_t analysisMode=mLocal)
 104 {
 105  // analysisMode: if analysisMode = mLocal -> analyze data on your computer using aliroot
 106  //               if analysisMode = mLocalSource -> analyze data on your computer using root + source files
 107
 108  // Cross-check user settings:
 109  CrossCheckSettings();
 110
 111  // Load needed libraries:
 112  LoadLibrariesPC(analysisMode);
 113
 114  // Global settings which will affect all plots:
 115  GlobalSettings();
 116
 117  // Access all common output files:
 118  AccessCommonOutputFiles(commonOutputFileName);
 119
 120  // Get from common output files the lists holding histograms for each method:
 121  GetLists();
 122
 123  // Get histograms with results for cumulants vs multiplicity:
 124  GetHistograms();
 125
 126  // Determine ranges for plots:
 127  DetermineMinMax();
 128
 129  // Make lines which indicate theoretical contributions of flow to cumulants:
 130  Lines();
 131
 132  // Print number of events and average multiplicities for each common output file:
 133  Print();
 134
 135  // Make plots for cumulants vs multiplicity: // to be improved
 136  PlotCumulantsVsM();
 137
 138  // Make plots for reference flow vs multiplicity:
 139  if(showRefFlowVsM){PlotRefFlowVsM();}
 140
 141 } // end of void plotCumulants(Int_t analysisMode=mLocal)
 142
 143 // =====================================================================================
 144
 145 void PlotRefFlowVsM()
 146 {
 147  // Make plots for reference flow vs multiplicity.
 148
 149  // Calculate reference flow from cumulants vs M:
 150  CalculateReferenceFlowFromCumulantsVsM();
 151
 152  TCanvas *cRefFlowVsM = new TCanvas("cRefFlowVsM","Reference Flow");
 153  cRefFlowVsM->Divide(2,2);
 154
 155  TLegend *lRefFlowVsM = new TLegend(0.1,0.7,0.33,0.9);
 156  lRefFlowVsM->SetFillStyle(0);
 157  lRefFlowVsM->SetHeader("     Therminator (20-30%)");
 158
 159  TString refFlowVsMFlag[4] = {"v_{2}{2,QC}","v_{2}{4,QC}","v_{2}{6,QC}","v_{2}{8,QC}"};
 160
 161  for(Int_t f=0;f<nFiles;f++)
 162  {
 163   for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++)
 164   {
 165    for(Int_t co=0;co<4;co++) // cumulant order
 166    {
 167     cRefFlowVsM->cd(co+1);
 168     // Style histogram:
 169     if(f==0) // superimpose histograms from other files on top of this one
 170     {
 171      TH1D *styleHist = (TH1D*)StyleHist(refFlowVsMFlag[co].Data(),co)->Clone();
 172      if(styleHist)
 173      {
 174       styleHist->GetXaxis()->SetRangeUser(refFlowVsMxRange[0],refFlowVsMxRange[1]);
 175       styleHist->GetYaxis()->SetRangeUser(refFlowVsMyRange[0],refFlowVsMyRange[1]);
 176       styleHist->Draw();
 177      }
 178     } // end of if(f==0)
 179     // Plot first the mesh for reference flow of specified order in all pads:
 180     if(refFlowVsMeshOrder != -1)
 181     {
 182      TGraph *rfMeshOutOf = GetErrorMesh(refFlowVsM[f][m][refFlowVsMeshOrder]);
 183      if(rfMeshOutOf)
 184      {
 185       rfMeshOutOf->SetFillColor(refFlowVsMMeshColor);
 186       rfMeshOutOf->Draw("lfsame");
 187       if(refFlowVsMeshOrder==co){lRefFlowVsM->AddEntry(rfMeshOutOf,Form("v_{2}{%i,QC} stat. error",2*(refFlowVsMeshOrder+1)),"f");}
 188      }
 189     } // end of if(refFlowVsMeshOrder != -1)
 190     // Plot reference flow vs M:
 191     refFlowVsM[f][m][co]->SetMarkerStyle(refFlowVsMMarkerStyle[f]);
 192     refFlowVsM[f][m][co]->SetMarkerColor(refFlowVsMMarkerColor[co]);
 193     refFlowVsM[f][m][co]->Draw("e1same");
 194     lRefFlowVsM->AddEntry(refFlowVsM[f][m][co],Form("v_{2}{%i,QC}, %s",2*(co+1),files[f].Data()),"p");
 195    } // end of for(Int_t co=0;co<4;co++) // cumulant order
 196   } // end of for(Int_t m=0;m<nMethods;m++)
 197  } // end of for(Int_t f=0;f<nFiles;f++)
 198
 199  // Draw a common legend in the 1st pad:
 200  cRefFlowVsM->cd(1);
 201  lRefFlowVsM->Draw("same");
 202
 203 } // end of void PlotRefFlowVsM()
 204
 205 // =====================================================================================
 206
 207 void PlotCumulantsVsM()
 208 {
 209  // Make plots for cumulants vs multiplicity: // to be improved
 210
 211  TCanvas *c = NULL;
 212  Int_t coMax = 0;
 213  if(!plotOnly2ndAnd4thOrderCumulant)
 214  {
 215   c = new TCanvas("c","cumulants");
 216   c->Divide(2,2);
 217   coMax = 4;
 218  } else
 219    {
 220     c = new TCanvas("c","cumulants",1200,500);
 221     c->Divide(2,1);
 222     coMax = 2;
 223    }
 224
 225  TLegend *legend = new TLegend(0.1,0.7,0.33,0.9);
 226  legend->SetFillStyle(0);
 227  //legend->SetHeader("     minClustersTpcRP");
 228
 229  TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
 230
 231  for(Int_t co=0;co<coMax;co++) // cumulant order
 232  {
 233   c->cd(co+1);
 234   TH1D *styleHist = (TH1D*) StyleHist(qcFlag[co].Data(),co)->Clone();
 235   if(styleHist){styleHist->Draw();}
 236   if(co==0)
 237   {
 238    if(showToyModel)
 239    {
 240     for(Int_t ntm=0;ntm<nToyModels;ntm++)
 241     {
 242      TF1 *tm = ToyModel(co,k[ntm],vn[ntm],v2n[ntm]);
 243      if(tm)
 244      {
 245       tm->SetLineColor(lineColorToyModel[ntm]);
 246       tm->Draw("same");
 247      }
 248     } // for(Int_t ntm=0;ntm<nToyModels;ntm++)
 249    } // end of if(showToyModel)
 250   } // end of if(co==0)
 251   if(co==1)
 252   {
 253    if(showToyModel)
 254    {
 255     for(Int_t ntm=0;ntm<nToyModels;ntm++)
 256     {
 257      TF1 *tm = ToyModel(co,k[ntm],vn[ntm],v2n[ntm]);
 258      if(tm)
 259      {
 260       tm->SetLineColor(lineColorToyModel[ntm]);
 261       tm->Draw("same");
 262      }
 263     } // for(Int_t ntm=0;ntm<nToyModels;ntm++)
 264    } // end of if(showToyModel)
 265   } // end of if(co==1)
 266
 267   // simulations:
 268   for(Int_t s=0;s<nSim;s++)
 269   {
 270    // Monte Carlo QC:
 271    if(showAlsoMCEPResults)
 272    {
 273     TGraph *mcepQC = GetErrorMesh(cumulantsVsM[s][2][co]);
 274     if(mcepQC)
 275     {
 276      mcepQC->SetFillColor(mcepMeshColor[s]);
 277      mcepQC->Draw("lfsame");
 278     }
 279     if(co==0){legend->AddEntry(mcepQC,Form("%s (MC)",legendEntry[s].Data()),"f");}
 280    } // end of if(showAlsoMCEPResults)
 281    // QC:
 282    TGraph *errorMesh = GetErrorMesh(cumulantsVsM[s][0][co]);
 283    if(errorMesh)
 284    {
 285     errorMesh->SetFillColor(meshColor[s]);
 286     if(!(showOnlyMCEPResults && showAlsoMCEPResults)){errorMesh->Draw("lfsame");}
 287    }
 288    if(co==0 && !(showOnlyMCEPResults && showAlsoMCEPResults)){legend->AddEntry(errorMesh,legendEntry[s].Data(),"f");}
 289   } // end of if(Int_t s=0;s<nSim;s++)
 290   // Theoretical lines:
 291   if(showTheoreticalLines)
 292   {
 293    for(Int_t fv=0;fv<nFlowValues;fv++)
 294    {
 295     lines[fv][co]->Draw("lsame");
 296     if(co==0){legend->AddEntry(lines[fv][co],Form("v_{2} = %g",v[fv]),"l");}
 297    }
 298   }
 299   // data:
 300   for(Int_t f=nSim;f<nFiles;f++)
 301   {
 302    // QC results:
 303    if(cumulantsVsM[f][0][co])
 304    {
 305     cumulantsVsM[f][0][co]->Draw("e1same");
 306     cumulantsVsM[f][0][co]->SetMarkerStyle(markerStyle[f-nSim]);
 307     cumulantsVsM[f][0][co]->SetMarkerColor(markerColor[f-nSim]);
 308     if(co==0)
 309     {
 310      if(showAlsoGFCResults)
 311      {
 312       legend->AddEntry(cumulantsVsM[f][0][co],Form("%s (QC)",legendEntry[f].Data()),"p");
 313      } else
 314        {
 315         legend->AddEntry(cumulantsVsM[f][0][co],legendEntry[f].Data(),"p");
 316        }
 317     }
 318    }
 319    // GFC results:
 320    if(showAlsoGFCResults && cumulantsVsM[f][1][co])
 321    {
 322     cumulantsVsM[f][1][co]->Draw("lsame");
 323     cumulantsVsM[f][1][co]->SetLineStyle(gfcLineStyle);
 324     if(co==0){legend->AddEntry(cumulantsVsM[f][1][co],Form("%s (GFC)",legendEntry[f].Data()),"l");}
 325    }
 326   }
 327   // Draw legend:
 328   if(co==0){legend->Draw("same");}
 329  } // end of for(Int_t co=0;co<4;co++) // cumulant order
 330
 331  // Plot also <reference multiplicity> vs # of RPs:
 332  if(plotCumulantsVsReferenceMultiplicity && showReferenceMultiplicityVsNoOfRPs)
 333  {
 334   TCanvas *cRefMultVsNoOfRPs = new TCanvas("cRefMultVsNoOfRPs","#LTreference multiplicity#GT vs # of RPs",1200,600);
 335   cRefMultVsNoOfRPs->Divide(nFiles,1);
 336   for(Int_t f=0;f<nFiles;f++)
 337   {
 338    cRefMultVsNoOfRPs->cd(f+1);
 339    if(refMultVsNoOfRPs[f])
 340    {
 341     refMultVsNoOfRPs[f]->SetTitle(legendEntry[f].Data());
 342     refMultVsNoOfRPs[f]->GetXaxis()->SetRangeUser(0,refMultVsNoOfRPs[f]->FindLastBinAbove());
 343     refMultVsNoOfRPs[f]->Draw();
 344    } else
 345      {
 346       cout<<endl;
 347       cout<<"WARNING: refMultVsNoOfRPs[f] is NULL in Plot(), f = "<<f<<" !!!!"<<endl;
 348       cout<<endl;
 349      }
 350   }
 351  }
 352
 353 } // end of void PlotCumulantsVsM()
 354
 355 // =====================================================================================
 356
 357 void CrossCheckSettings()
 358 {
 359  // Cross-check user settings in this method.
 360
 361  if(showAlsoGFCResults && rebin)
 362  {
 363   cout<<endl;
 364   cout<<" WARNING: Rebinning in M not supported for GFC yet !!!!"<<endl;
 365   cout<<endl;
 366   exit(0);
 367  }
 368
 369  if(showAlsoGFCResults && plotCumulantsVsReferenceMultiplicity)
 370  {
 371   cout<<endl;
 372   cout<<" WARNING: Showing GFC versus <reference multiplicity> not supported yet !!!!"<<endl;
 373   cout<<endl;
 374   exit(0);
 375  }
 376
 377 } // end of void CrossCheckSettings()
 378
 379 // =====================================================================================
 380
 381 void Lines()
 382 {
 383  // Make lines denoting theoretical contribution of flow to cumulants.
 384
 385  for(Int_t fv=0;fv<nFlowValues;fv++)
 386  {
 387   lines[fv][0] = new TGraph(2);
 388   lines[fv][0]->SetPoint(0,xMin[0],pow(v[fv],2));
 389   lines[fv][0]->SetPoint(1,xMax[0]+0.5,pow(v[fv],2));
 390   lines[fv][0]->SetLineColor(lineColor[fv]);
 391   lines[fv][1] = new TGraph(2);
 392   lines[fv][1]->SetPoint(0,xMin[1],-pow(v[fv],4));
 393   lines[fv][1]->SetPoint(1,xMax[1]+0.5,-pow(v[fv],4));
 394   lines[fv][1]->SetLineColor(lineColor[fv]);
 395   lines[fv][2] = new TGraph(2);
 396   lines[fv][2]->SetPoint(0,xMin[2],4.*pow(v[fv],6));
 397   lines[fv][2]->SetPoint(1,xMax[2]+0.5,4.*pow(v[fv],6));
 398   lines[fv][2]->SetLineColor(lineColor[fv]);
 399   lines[fv][3] = new TGraph(2);
 400   lines[fv][3]->SetPoint(0,xMin[3],-33.*pow(v[fv],8));
 401   lines[fv][3]->SetPoint(1,xMax[3]+0.5,-33.*pow(v[fv],8));
 402   lines[fv][3]->SetLineColor(lineColor[fv]);
 403  }
 404
 405 } // end of void Lines()
 406
 407 // =====================================================================================
 408
 409 TF1* ToyModel(Int_t co, Double_t k, Double_t vn, Double_t v2n)
 410 {
 411  // Make theoretical curves for the toy model.
 412
 413  TF1 *function = NULL;
 414
 415  if(co==0)
 416  {
 417   function = new TF1("function","([1]*[1]*(x-0.5-[0])+[0]-1)/((x-0.5)-1)",1.5,1000);
 418   function->SetParameter(0,k);
 419   function->SetParameter(1,vn);
 420  }
 421  else if(co==1)
 422  {
 423   function = new TF1("function","-1.*([1]*[1]*[1]*[1]*([0]-(x-0.5))*(12.*[0]-6.*[0]*[0]-12.*(x-0.5)-5.*[0]*(x-0.5)+6.*[0]*[0]*(x-0.5)+9.*(x-0.5)*(x-0.5)-3.*[0]*(x-0.5)*(x-0.5)-(x-0.5)*(x-0.5)*(x-0.5))+[1]*[1]*4.*([0]-1.)*(x-0.5-[0])*([0]*(x-0.5-1.)-2.*(x-0.5-2.))-[1]*[1]*[2]*2.*(x-0.5-1.)*([0]-1.)*(x-0.5-[0])*(x-0.5-2.*[0])-[2]*[2]*([0]-1.)*([0]-1.)*((x-0.5)-[0])*(x-0.5-1.)+([0]-1.)*(-6.+9*[0]-[0]*[0]+2.*(x-0.5)-5.*[0]*(x-0.5)+[0]*[0]*(x-0.5)))/((x-0.5-1.)*(x-0.5-1.)*(x-0.5-2.)*(x-0.5-3.))",3.5,100);
 424   function->SetParameter(0,k);
 425   function->SetParameter(1,vn);
 426   function->SetParameter(2,v2n);
 427  }
 428
 429  return function;
 430
 431 } // end of TF1* ToyModel(Int_t co)
 432
 433 // =====================================================================================
 434
 435 void CalculateReferenceFlowFromCumulantsVsM()
 436 {
 437  // Calculate reference flow from cumulants vs M:
 438
 439  for(Int_t f=0;f<nFiles;f++)
 440  {
 441   for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++)
 442   {
 443    for(Int_t co=0;co<4;co++) // cumulant order
 444    {
 445     if(cumulantsVsM[f][m][co])
 446     {
 447      refFlowVsM[f][m][co] = (TH1D*)cumulantsVsM[f][m][co]->Clone(Form("%i,%i,%i",f,m,co));
 448      if(!refFlowVsM[f][m][co]){cout<<" WARNING: "<<Form("refFlowVsM[%i][%i][%i]",f,m,co)<<" is NULL in PlotRefFlowVsM() !!!!"<<exit(0)<<endl;}
 449      refFlowVsM[f][m][co]->Reset(); // to have empty histogram with the same binning as the one with cumulants
 450      Int_t nBins = refFlowVsM[f][m][co]->GetNbinsX();
 451      for(Int_t b=1;b<=nBins;b++)
 452      {
 453       Double_t qcVsM = cumulantsVsM[f][m][co]->GetBinContent(b); // QC vs M
 454       Double_t qcErrorVsM = cumulantsVsM[f][m][co]->GetBinError(b); // QC vs M stat. error
 455       Double_t vVsM = 0.; // reference flow vs M
 456       Double_t vErrorVsM = 0.; // reference flow vs M stat. error
 457       if(co==0) // 2nd order
 458       {
 459        if(qcVsM>0.)
 460        {
 461         vVsM = pow(qcVsM,1./2.);
 462         vErrorVsM = (1./2.)*pow(qcVsM,-1./2.)*qcErrorVsM;
 463        }
 464       } // end of if(co==0) 2nd order
 465       else if(co==1) // 4th order
 466       {
 467        if(qcVsM<0.)
 468        {
 469         vVsM = pow(-1.*qcVsM,1./4.);
 470         vErrorVsM = (1./4.)*pow(-qcVsM,-3./4.)*qcErrorVsM;
 471        }
 472       } // end of if(co==1) 4th order
 473       else if(co==2) // 6th order
 474       {
 475        if(qcVsM>0.)
 476        {
 477         vVsM = pow((1./4.)*qcVsM,1./6.);
 478         vErrorVsM = (1./6.)*pow(2.,-1./3.)*pow(qcVsM,-5./6.)*qcErrorVsM;
 479        }
 480       } // end of if(co==2) 6th order
 481       else if(co==3) // 8th order
 482       {
 483        if(qcVsM<0.)
 484        {
 485         vVsM = pow((-1./33.)*qcVsM,1./8.);
 486         vErrorVsM = (1./8.)*pow(33.,-1./8.)*pow(-qcVsM,-7./8.)*qcErrorVsM;
 487        }
 488       } // end of if(co==3) 8th order
 489       // Store final results and statisticial errror for reference flow:
 490       refFlowVsM[f][m][co]->SetBinContent(b,vVsM);
 491       refFlowVsM[f][m][co]->SetBinError(b,vErrorVsM);
 492      } // end of for(Int_t b=1;b<=nBins;b++)
 493     } else
 494       {
 495        cout<<endl;
 496        cout<<" WARNING: "<<Form("cumulantsVsM[%i][%i][%i]",f,m,co)<<" is NULL in CalculateReferenceFlowFromCumulantsVsM() !!!!"<<endl;
 497        cout<<endl;
 498       }
 499    } // end of for(Int_t co=0;co<4;co++) // cumulant order
 500   } // end of for(Int_t m=0;m<nMethods;m++)
 501  } // end of for(Int_t f=0;f<nFiles;f++)
 502
 503 } // end of void CalculateReferenceFlowFromCumulantsVsM()
 504
 505 // =====================================================================================
 506
 507 void Print()
 508 {
 509  // Print number of events and average multiplicities for each common output file.
 510
 511  cout<<endl;
 512  cout<<"Accessed files:"<<endl;
 513  cout<<endl;
 514  for(Int_t f=0;f<nFiles;f++)
 515  {
 516   cout<<commonOutputFiles[f]->GetName()<<endl;
 517   for(Int_t m=0;m<nMethods;m++)
 518   {
 519    AliFlowCommonHist *commonHist = NULL;
 520    if(lists[f][m])
 521    {
 522     commonHist = dynamic_cast<AliFlowCommonHist*> (lists[f][m]->FindObject(Form("AliFlowCommonHist%s",method[m].Data())));
 523    }
 524    Double_t nEvts = 0.;
 525    Double_t AvM = 0.;
 526    if(commonHist && commonHist->GetHistMultRP())
 527    {
 528     nEvts = commonHist->GetHistMultRP()->GetEntries();
 529     AvM = commonHist->GetHistMultRP()->GetMean();
 530    }
 531    if(!(strcmp(method[m].Data(),"QC")))
 532    {
 533     cout<<Form("%s:",method[m].Data())<<"  <M> = "<<AvM<<", N = "<<nEvts<<endl;
 534    }
 535    if(!(strcmp(method[m].Data(),"GFC")) && showAlsoGFCResults)
 536    {
 537     cout<<Form("%s:",method[m].Data())<<"  <M> = "<<AvM<<", N = "<<nEvts<<endl;
 538    }
 539   } // end of for(Int_t m=0;m<nMethods;m++)
 540   cout<<endl;
 541  } // end of for(Int_t f=0;f<nFiles;f++)
 542
 543 } // end of void Print()
 544
 545 // =====================================================================================
 546
 547 void DetermineMinMax()
 548 {
 549  // Determine ranges for plots.
 550
 551  for(Int_t co=0;co<4;co++)
 552  {
 553   xMin[co] = 0.; yMin[co] = 44.;
 554   xMax[co] = -440000.; yMax[co] = -44.;
 555  }
 556
 557  Double_t tfc[nFlowValues][4] = {{0.}}; // theoretical flow contribution
 558  for(Int_t fv=0;fv<nFlowValues;fv++)
 559  {
 560   tfc[fv][0] = pow(v[fv],2);
 561   tfc[fv][1] = -pow(v[fv],4);
 562   tfc[fv][2] = 4.*pow(v[fv],6);
 563   tfc[fv][3] = -33.*pow(v[fv],8);
 564  }
 565
 566  for(Int_t f=0;f<nFiles;f++)
 567  {
 568   for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++)
 569   {
 570    for(Int_t co=0;co<4;co++)
 571    {
 572     if(cumulantsVsM[f][m][co])
 573     {
 574      Int_t nBins = cumulantsVsM[f][m][co]->GetXaxis()->GetNbins();
 575      for(Int_t b=1;b<=nBins;b++)
 576      {
 577       Double_t result = cumulantsVsM[f][m][co]->GetBinContent(b);
 578       Double_t error = cumulantsVsM[f][m][co]->GetBinError(b);
 579       if(TMath::Abs(result)>1.e-44 && TMath::Abs(error)>1.e-44)
 580       {
 581        // y-axis:
 582        if(yMin[co] > result-error){yMin[co] = result-error;} // min value
 583        if(yMax[co] < result+error){yMax[co] = result+error;} // max value
 584        // x-axis:
 585        if(xMax[co] < cumulantsVsM[f][m][co]->GetBinLowEdge(b+1)){xMax[co] = cumulantsVsM[f][m][co]->GetBinLowEdge(b+1);}
 586       }
 587      } // end of for(Int_t b=1;b<=cumulantsVsM[f][m][co]->GetXaxis()->GetNbins();b++)
 588      // theoretical contributions:
 589      if(showTheoreticalLines)
 590      {
 591       for(Int_t fv=0;fv<nFlowValues;fv++)
 592       {
 593        if(yMin[co] > tfc[fv][co]) {yMin[co] = tfc[fv][co];} // min value
 594        if(yMax[co] < tfc[fv][co]) {yMax[co] = tfc[fv][co];} // max value
 595       } // end of for(Int_t fv=0;fv<nFlowValues;fv++)
 596      } // end of if(showTheoreticalLines)
 597     } // end of if(cumulantsVsM[f][m][co])
 598    } // end of for(Int_t co=0;co<4;co++)
 599   } // end of for(Int_t m=0;m<nMethods;m++)
 600  } // end of for(Int_t f=0;f<nFiles;f++)
 601
 602 } // end of void DetermineMinMax()
 603
 604 // =====================================================================================
 605
 606 void GetHistograms()
 607 {
 608  // Get all histograms and profiles.
 609
 610  // a) Get profiles holding results for <refMult> vs number of Reference Particles (RPs);
 611  // b) Get histograms holding results for cumulants and reference flow vs multiplicity.
 612
 613  // a) Get profiles holding results for <refMult> vs number of Reference Particles (RPs):
 614  if(plotCumulantsVsReferenceMultiplicity){GetProfileRefMultVsNoOfRPs();}
 615
 616  // b) Get histograms holding results for cumulants and reference flow vs multiplicity:
 617  TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
 618  TString gfcFlag[4] = {"GFC{2}","GFC{4}","GFC{6}","GFC{8}"};
 619  for(Int_t f=0;f<nFiles;f++)
 620  {
 621   for(Int_t m=0;m<nMethods;m++)
 622   {
 623    TList *temp = NULL;
 624    if(!(strcmp(method[m].Data(),"QC")) && lists[f][m])
 625    {
 626     temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Integrated Flow"));
 627     if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));}
 628     if(temp)
 629     {
 630      for(Int_t co=0;co<4;co++)
 631      {
 632       // Cumulants vs multiplicity:
 633       cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fIntFlowQcumulantsVsM, %s",qcFlag[co].Data())));
 634       if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co])
 635       {
 636        cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f);
 637       }
 638       if(rebin && cumulantsVsM[f][m][co])
 639       {
 640        cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co);
 641       }
 642      } // end of for(Int_t co=0;co<4;co++)
 643     } // end of if(temp)
 644    } // end of if(!(strcmp(method[m].Data(),"QC")))
 645    else if(!(strcmp(method[m].Data(),"GFC")) && lists[f][m])
 646    {
 647     temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Reference Flow"));
 648     if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));}
 649     if(temp)
 650     {
 651      for(Int_t co=0;co<4;co++)
 652      {
 653       cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fReferenceFlowCumulantsVsM, %s",gfcFlag[co].Data())));
 654       if(showAlsoGFCResults && !cumulantsVsM[f][m][co])
 655       {
 656        cout<<endl;
 657        cout<<Form(" WARNING: Couldn't access histogram fReferenceFlowCumulantsVsM, %s in the ",gfcFlag[co].Data())<<endl;
 658        cout<<"          file "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
 659        cout<<"          Did you enable calculation of GFC vs M in the analysis which produced this file?"<<endl;
 660        cout<<endl;
 661       }
 662       if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co])
 663       {
 664        cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f);
 665       }
 666       if(rebin && cumulantsVsM[f][m][co])
 667       {
 668        cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co);
 669       }
 670      } // end of for(Int_t co=0;co<4;co++)
 671     } // end of if(temp)
 672    } // end of else if(!(strcmp(method[m].Data(),"GFC")))
 673    else if(!(strcmp(method[m].Data(),"MCEP")) && lists[f][m])
 674    {
 675     TProfile *mcepVsM = dynamic_cast<TProfile*> (lists[f][m]->FindObject("FlowPro_VsM_MCEP"));
 676     if(!mcepVsM)
 677     {
 678      cout<<endl;
 679      cout<<" WARNING: Couldn't access TProfile FlowPro_VsM_MCEP in the file "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
 680      cout<<endl;
 681      return;
 682     }
 683     for(Int_t co=0;co<4;co++)
 684     {
 685      cumulantsVsM[f][m][co] = new TH1D("",Form("Monte Carlo QC{%i} #font[72]{vs} M",2*(co+1)),
 686                                        mcepVsM->GetNbinsX(),mcepVsM->GetBinLowEdge(1),mcepVsM->GetBinLowEdge(1+mcepVsM->GetNbinsX()));
 687     } // end of for(Int_t co=0;co<4;co++)
 688     for(Int_t b=1;b<=mcepVsM->GetNbinsX();b++)
 689     {
 690      Double_t v = mcepVsM->GetBinContent(b);
 691      Double_t vError = mcepVsM->GetBinError(b);
 692      if(TMath::Abs(v)<1.e-44 || TMath::Abs(vError)<1.e-44){continue;}
 693      // Monte Carlo QC{2}:
 694      Double_t qc2 = pow(v,2.);
 695      Double_t qc2Error = 2.*TMath::Abs(v)*vError; // error propagation for f(x) = x^2
 696      cumulantsVsM[f][m][0]->SetBinContent(b,qc2);
 697      cumulantsVsM[f][m][0]->SetBinError(b,qc2Error);
 698      // Monte Carlo QC{4}:
 699      Double_t qc4 = -pow(v,4.);
 700      Double_t qc4Error = 4.*TMath::Abs(pow(v,3.))*vError; // error propagation for f(x) = -x^4
 701      cumulantsVsM[f][m][1]->SetBinContent(b,qc4);
 702      cumulantsVsM[f][m][1]->SetBinError(b,qc4Error);
 703      // Monte Carlo QC{6}:
 704      Double_t qc6 = 4.*pow(v,6.);
 705      Double_t qc6Error = 24.*TMath::Abs(pow(v,5.))*vError; // error propagation for f(x) = 4x^6
 706      cumulantsVsM[f][m][2]->SetBinContent(b,qc6);
 707      cumulantsVsM[f][m][2]->SetBinError(b,qc6Error);
 708      // Monte Carlo QC{8}:
 709      Double_t qc8 = -33.*pow(v,8.);
 710      Double_t qc8Error = 264.*TMath::Abs(pow(v,7.))*vError; // error propagation for f(x) = -33x^8
 711      cumulantsVsM[f][m][3]->SetBinContent(b,qc8);
 712      cumulantsVsM[f][m][3]->SetBinError(b,qc8Error);
 713     } // end of for(Int_t b=1;b<=mcepVsM->GetNbinsX();b++)
 714     for(Int_t co=0;co<4;co++)
 715     {
 716      if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co])
 717      {
 718       cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f);
 719      }
 720      if(rebin && cumulantsVsM[f][m][co])
 721      {
 722       cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co);
 723      }
 724     } // end of for(Int_t co=0;co<4;co++)
 725    } // end of else if(!(strcmp(method[m].Data(),"MCEP")))
 726   } // end of  for(Int_t m=0;m<nMethods;m++)
 727  } // end of for(Int_t f=0;f<nFiles;f++)
 728
 729  Int_t counter = 0;
 730  for(Int_t f=0;f<nFiles;f++)
 731  {
 732   for(Int_t m=0;m<nMethods;m++)
 733   {
 734    for(Int_t co=0;co<4;co++)
 735    {
 736     if(cumulantsVsM[f][m][co]){counter++;}
 737    }
 738   }
 739  }
 740
 741  if(counter == 0)
 742  {
 743   cout<<endl;
 744   cout<<" WARNING: Couldn't access a single histogram with results vs multiplicity !!!!"<<endl;
 745   cout<<"          Did you enable this calculation before running over data?"<<endl;
 746   cout<<endl;
 747   exit(0);
 748  }
 749
 750 } // end of void GetHistograms()
 751
 752 // =====================================================================================
 753
 754 void GetProfileRefMultVsNoOfRPs()
 755 {
 756  // Get profiles holding results for <reference multiplicity> vs number of Reference Particles (RPs).
 757
 758  // Set here from which method's output file this profile will be accessed:
 759  TString methodName = "QC"; // Alternatives are GFC and MCEP
 760  Int_t i = -1;
 761  for(Int_t m=0;m<nMethods;m++)
 762  {
 763   if(method[m] == methodName){i=m;}
 764  }
 765  if(i==-1)
 766  {
 767   cout<<endl;
 768   cout<<" WARNING: Unknown method name in GetProfileRefMultVsNoOfRPs() !!!!"<<endl;
 769   cout<<"          Try something else for TString methodName in GetProfileRefMultVsNoOfRPs()."<<endl;
 770   cout<<endl;exit(0);
 771  }
 772
 773  for(Int_t f=0;f<nFiles;f++)
 774  {
 775   AliFlowCommonHist *commonHist = NULL;
 776   if(lists[f][i])
 777   {
 778    commonHist = dynamic_cast<AliFlowCommonHist*> (lists[f][i]->FindObject(Form("AliFlowCommonHist%s",method[i].Data())));
 779   } else
 780     {
 781      cout<<endl;
 782      cout<<Form(" WARNING: lists[%i][%i] is NULL in GetProfileRefMultVsNoOfRPs() !!!!",f,i)<<endl;
 783      cout<<Form("          Did you use method %s in the analysis which has produced output",method[i].Data())<<endl;
 784      cout<<Form("          file %s?",commonOutputFiles[f]->GetName())<<endl;
 785     }
 786   if(commonHist && commonHist->GetRefMultVsNoOfRPs())
 787   {
 788    refMultVsNoOfRPs[f] = commonHist->GetRefMultVsNoOfRPs();
 789   } else
 790     {
 791      cout<<endl;
 792      cout<<" WARNING: commonHist && commonHist->GetRefMultVsNoOfRPs() is NULL in GetProfileRefMultVsNoOfRPs() !!!!"<<endl;
 793      cout<<endl;
 794     }
 795  } // end of for(Int_t f=0;f<nFiles;f++)
 796
 797  return;
 798
 799 } // end of void GetProfileRefMultVsNoOfRPs()
 800
 801 // =====================================================================================
 802
 803 TH1D* Map(TH1D *hist, Int_t f)
 804 {
 805  // Map cumulant versus # of RPs into cumulants versus <reference multiplicity>.
 806
 807  if(!refMultVsNoOfRPs[f])
 808  {
 809   cout<<endl;
 810   cout<<Form(" WARNING: refMultVsNoOfRPs[%i] is NULL in Map(...) !!!!",f)<<endl;
 811   cout<<endl;
 812  }
 813
 814  Int_t rpMinBin = refMultVsNoOfRPs[f]->FindFirstBinAbove(); // FindFirstBinAbove(Double_t threshold = 0, Int_t axis = 1)
 815  Int_t rpMaxBin = refMultVsNoOfRPs[f]->FindLastBinAbove(); // FindLastBinAbove(Double_t threshold = 0, Int_t axis = 1)
 816  Int_t rmMinBin = 440000;;
 817  Int_t rmMaxBin = (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetMaximum());
 818  for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // non-empty # of RPs bins
 819  {
 820   if(refMultVsNoOfRPs[f]->GetBinContent(rpBin)>0. && refMultVsNoOfRPs[f]->GetBinContent(rpBin)<rmMinBin)
 821   {
 822    rmMinBin = (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin));
 823   }
 824  } // end of for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // non-empty # of RPs bins
 825
 826  if(hist)
 827  {
 828   temp = (TH1D*) hist->Clone();
 829   temp->Reset();
 830   for(Int_t rmBin=rmMinBin;rmBin<=rmMaxBin;rmBin++) // reference multiplicity bins
 831   {
 832    Double_t value = 0.;
 833    Double_t error = 0.;
 834    Double_t dSum1 = 0.; // sum value_i/(error_i)^2
 835    Double_t dSum2 = 0.; // sum 1/(error_i)^2
 836    for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // # of RPs bins
 837    {
 838     if((Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin)) >= temp->GetBinLowEdge(rmBin) &&
 839        (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin)) < temp->GetBinLowEdge(rmBin+1))
 840     {
 841      value = hist->GetBinContent(rpBin);
 842      error = hist->GetBinError(rpBin);
 843      if(error>0.)
 844      {
 845       dSum1+=value/(error*error);
 846       dSum2+=1./(error*error);
 847      }
 848     }
 849    } // end of for(Int_t rpBin=1;rpBin<=nBins;rpBin++) // # of RPs bins
 850    if(dSum2>0.)
 851    {
 852     temp->SetBinContent(rmBin,dSum1/dSum2);
 853     temp->SetBinError(rmBin,pow(1./dSum2,0.5));
 854    }
 855   } // end of for(Int_t rmBin=1;rmBin<=nBins;rmBin++) // reference multiplicity bins
 856  } // end of if(hist)
 857
 858  return temp;
 859
 860 } // end of Map()
 861
 862 // =====================================================================================
 863
 864 TH1D* Rebin(TH1D *hist, Int_t co)
 865 {
 866  // Rebin original histograms.
 867
 868  if(nMergedBins[co] == 0)
 869  {
 870   cout<<endl;
 871   cout<<Form(" WARNING: nMergedBins[%i] == 0 !!!!",co)<<endl;
 872   cout<<endl;
 873   exit(0);
 874  }
 875  if(!hist)
 876  {
 877   cout<<endl;
 878   cout<<" WARNING: hist is NULL in Rebin() !!!!"<<endl;
 879   cout<<endl;
 880   exit(0);
 881  }
 882
 883  Int_t nBinsOld = hist->GetXaxis()->GetNbins();
 884  if(nBinsOld == 0){cout<<" WARNING: nBinsOld == 0 !!!!"<<endl;exit(0);}
 885  Double_t xMinOld = hist->GetXaxis()->GetXmin();
 886  Double_t xMaxOld = hist->GetXaxis()->GetXmax();
 887  Double_t binWidthOld = (xMaxOld-xMinOld)/nBinsOld;
 888  Int_t nBinsNew = TMath::Floor(nBinsOld/nMergedBins[co]);
 889  Double_t xMinNew = xMinOld;
 890  Double_t xMaxNew = xMinOld + nBinsNew*nMergedBins[co]*binWidthOld;
 891
 892  TH1D *temp = new TH1D("","",nBinsNew,xMinNew,xMaxNew); // rebinned histogram
 893  Int_t binNew = 1;
 894  Double_t value = 0.;
 895  Double_t error = 0.;
 896  Double_t dSum1 = 0.; // sum value_i/(error_i)^2
 897  Double_t dSum2 = 0.; // sum 1/(error_i)^2
 898  for(Int_t b=1;b<=nBinsOld;b++)
 899  {
 900   value = hist->GetBinContent(b);
 901   error = hist->GetBinError(b);
 902   if(error>0.)
 903   {
 904    dSum1+=value/(error*error);
 905    dSum2+=1./(error*error);
 906   }
 907   if(b%nMergedBins[co] == 0)
 908   {
 909    if(dSum2>0.)
 910    {
 911     temp->SetBinContent(binNew,dSum1/dSum2);
 912     temp->SetBinError(binNew,pow(1./dSum2,0.5));
 913    }
 914    binNew++;
 915    dSum1 = 0.;
 916    dSum2 = 0.;
 917   } // end of if(b%nMergedBins[co] == 0)
 918  } // end of for(Int_t b=1;b<=nBinsOld;b++)
 919
 920  return temp;
 921
 922 } // end of TH1D* Rebin(TH1D *hist, Int_t co)
 923
 924 // =====================================================================================
 925
 926 TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc)
 927 {
 928  TGraphErrors *ge = new TGraphErrors(nFiles);
 929  for(Int_t f=0;f<nFiles;f++)
 930  {
 931   ge->SetPoint(f,f+0.5,qc[f]->GetBinContent(bin+1));
 932   ge->SetPointError(f,0,qc[f]->GetBinError(bin+1));
 933  }
 934
 935  return ge;
 936
 937 } // end of TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc)
 938
 939 // =====================================================================================
 940
 941 TH1D* StyleHist(TString yAxisTitle, Int_t co)
 942 {
 943  // Style histogram.
 944
 945  TH1D *styleHist = new TH1D(yAxisTitle.Data(),"",(Int_t)xMax[co],0,xMax[co]);
 946  // y-axis:
 947  styleHist->GetYaxis()->SetRangeUser(yMin[co],yMax[co]);
 948
 949  if(plotCumulantsVsReferenceMultiplicity)
 950  {
 951   styleHist->GetXaxis()->SetTitle("#LTreference multiplicity#GT");
 952  } else
 953    {
 954     styleHist->GetXaxis()->SetTitle("# of RPs");
 955    }
 956
 957  styleHist->GetYaxis()->SetTitle(yAxisTitle.Data());
 958
 959  return styleHist;
 960
 961 } // end of TH1D* StyleHist(TString yAxisTitle, Int_t co)
 962
 963 // ===========================================================================================
 964
 965 TGraph* GetErrorMesh(TH1D *hist)
 966 {
 967  // Error mesh.
 968
 969  if(hist)
 970  {
 971   Int_t nBins = hist->GetNbinsX();
 972   Double_t value = 0.;
 973   Double_t error = 0.;
 974   // Counting the non-empty bins:
 975   Int_t nNonEmptyBins = 0;
 976   for(Int_t i=1;i<nBins+1;i++)
 977   {
 978    value = hist->GetBinContent(i);
 979    error = hist->GetBinError(i);
 980    if(TMath::Abs(value)>0.0 && error>0.0)
 981    {
 982     nNonEmptyBins++;
 983    }
 984   } // end of for(Int_t i=1;i<nBins+1;i++)
 985   // Error mesh:
 986   TGraph *errorMesh = new TGraph(2*nNonEmptyBins+1);
 987   Int_t count = 0;
 988   Double_t binCenter = 0.;
 989   for(Int_t i=1;i<nBins+1;i++)
 990   {
 991    value = hist->GetBinContent(i);
 992    error = hist->GetBinError(i);
 993    // Setting up the the mesh:
 994    if(TMath::Abs(value)>0.0 && error>0.0)
 995    {
 996     binCenter = hist->GetBinCenter(i);
 997     errorMesh->SetPoint(count,binCenter,value+error);
 998     errorMesh->SetPoint(2*nNonEmptyBins-(count+1),binCenter,value-error);
 999     count++;
1000    }
1001   } // end of for(Int_t i=1;i<nBins+1;i++)
1002   // Closing the mesh area:
1003   Double_t xStart = 0.;
1004   Double_t yStart = 0.;
1005   errorMesh->GetPoint(0,xStart,yStart);
1006   errorMesh->SetPoint(2*nNonEmptyBins,xStart,yStart);
1007  } // end if(hist)
1008
1009  return errorMesh;
1010
1011 } // end of TGraph* GetErrorMesh(TH1D *hist)
1012
1013 // ===========================================================================================
1014
1015 void GlobalSettings()
1016 {
1017  // Settings which will affect all plots.
1018
1019  gROOT->SetStyle("Plain"); // default color is white instead of gray
1020  gStyle->SetOptStat(0); // remove stat. box from all histos
1021  TGaxis::SetMaxDigits(4); // prefer exp notation for 5 and more significant digits
1022
1023 } // end of void GlobalSettings()
1024
1025 // ===========================================================================================
1026
1027 void GetLists()
1028 {
1029  // Get from common output files the lists holding histograms for each method.
1030
1031  TString fileName[nFiles][nMethods];
1032  TDirectoryFile *dirFile[nFiles][nMethods] = {{NULL}};
1033  TString listName[nFiles][nMethods];
1034  for(Int_t f=0;f<nFiles;f++)
1035  {
1036   for(Int_t i=0;i<nMethods;i++)
1037   {
1038    // Form a file name for each method:
1039    fileName[f][i]+="output";
1040    fileName[f][i]+=method[i].Data();
1041    fileName[f][i]+="analysis";
1042    fileName[f][i]+=type[f].Data();
1043    // Access this file:
1044    if(commonOutputFiles[f]){dirFile[f][i] = (TDirectoryFile*)commonOutputFiles[f]->FindObjectAny(fileName[f][i].Data());}
1045    // Form a list name for each method:
1046    if(dirFile[f][i])
1047    {
1048     TList* listTemp = dirFile[f][i]->GetListOfKeys();
1049     if(listTemp && listTemp->GetEntries() == 1)
1050     {
1051      listName[f][i] = listTemp->At(0)->GetName(); // to be improved - implemented better
1052      dirFile[f][i]->GetObject(listName[f][i].Data(),lists[f][i]);
1053     } else
1054       {
1055        cout<<endl;
1056        cout<<" WARNING: Couldn't find a list "<<listName[f][i].Data()<<" in "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
1057        cout<<"          Did you use method "<<method[i].Data()<<" in the analysis which produced this file?"<<endl;
1058       }
1059    } else
1060      {
1061       cout<<" WARNING: Couldn't find a file "<<fileName[f][i].Data()<<".root in "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl;
1062      }
1063   } // end of for(Int_t i=0;i<nMethods;i++)
1064  } // end of for(Int_t f=0;f<nFiles;f++)
1065
1066 } // end of void GetLists()
1067
1068 // ===========================================================================================
1069
1070 void AccessCommonOutputFiles(TString commonOutputFileName)
1071 {
1072  // Access all output files.
1073
1074  for(Int_t f=0;f<nFiles;f++)
1075  {
1076   if(bLookInSubdirectories)
1077   {
1078    if(!(gSystem->AccessPathName(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),kFileExists)))
1079    {
1080     commonOutputFiles[f] = TFile::Open(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),"READ");
1081    } else
1082      {
1083       cout<<endl;
1084       cout<<" WARNING: Couldn't find the file "<<Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data())<<" !!!!"<<endl;
1085       cout<<"          Did you specify correctly all paths in TString files[nFiles]?"<<endl;
1086       cout<<"          Or you are misusing 'Bool_t bLookInSubdirectories'? "<<endl;
1087       cout<<endl;
1088       exit(0);
1089      }
1090   } else // to if(bLookInSubdirectories)
1091     {
1092      if(!(gSystem->AccessPathName(Form("%s/%s%s",gSystem->pwd(),files[f].Data(),".root"),kFileExists)))
1093      {
1094       commonOutputFiles[f] = TFile::Open(Form("%s/%s%s",gSystem->pwd(),files[f].Data(),".root"),"READ");
1095      } else
1096        {
1097         cout<<endl;
1098         cout<<" WARNING: Couldn't find the file "<<Form("%s/%s%s",gSystem->pwd(),files[f].Data(),".root")<<" !!!!"<<endl;
1099         cout<<"          Did you specify correctly all file names in TString files[nFiles]?"<<endl;
1100         cout<<endl;
1101         exit(0);
1102        }
1103     }
1104  } // end of for(Int_t f=0;f<nFiles;f++)
1105
1106 } // void AccessCommonOutputFiles(TString commonOutputFileName);
1107
1108 // ===========================================================================================
1109
1110 void LoadLibrariesPC(const libModes analysisMode) {
1111
1112   //--------------------------------------
1113   // Load the needed libraries most of them already loaded by aliroot
1114   //--------------------------------------
1115   //gSystem->Load("libTree");
1116   gSystem->Load("libGeom");
1117   gSystem->Load("libVMC");
1118   gSystem->Load("libXMLIO");
1119   gSystem->Load("libPhysics");
1120
1121   //----------------------------------------------------------
1122   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
1123   //----------------------------------------------------------
1124   if (analysisMode==mLocal) {
1125     //--------------------------------------------------------
1126     // If you want to use already compiled libraries
1127     // in the aliroot distribution
1128     //--------------------------------------------------------
1129
1130     //==================================================================================
1131     //load needed libraries:
1132     gSystem->AddIncludePath("-I$ROOTSYS/include");
1133     //gSystem->Load("libTree");
1134
1135     // for AliRoot
1136     gSystem->AddIncludePath("-I$ALICE_ROOT/include");
1137     //gSystem->Load("libANALYSIS");
1138     gSystem->Load("libPWGflowBase");
1139     //cerr<<"libPWGflowBase loaded ..."<<endl;
1140
1141   }
1142
1143   else if (analysisMode==mLocalSource) {
1144
1145     // In root inline compile
1146
1147     // Constants
1148     gROOT->LoadMacro("Base/AliFlowCommonConstants.cxx+");
1149     gROOT->LoadMacro("Base/AliFlowLYZConstants.cxx+");
1150
1151     // Flow event
1152     gROOT->LoadMacro("Base/AliFlowVector.cxx+");
1153     gROOT->LoadMacro("Base/AliFlowTrackSimple.cxx+");
1154     gROOT->LoadMacro("Base/AliFlowTrackSimpleCuts.cxx+");
1155     gROOT->LoadMacro("Base/AliFlowEventSimple.cxx+");
1156
1157     // Output histosgrams
1158     gROOT->LoadMacro("Base/AliFlowCommonHist.cxx+");
1159     gROOT->LoadMacro("Base/AliFlowCommonHistResults.cxx+");
1160     gROOT->LoadMacro("Base/AliFlowLYZHist1.cxx+");
1161     gROOT->LoadMacro("Base/AliFlowLYZHist2.cxx+");
1162
1163     cout << "finished loading macros!" << endl;
1164
1165   }
1166
1167 } // end of void LoadLibraries(const libModes analysisMode)