PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /*************************************************
  19  * Flow analysis with cumulants. In this class   *
  20  * cumulants are calculated by making use of the *
  21  * formalism of generating functions proposed by *
  22  * Ollitrault et al.                             *
  23  *                                               *
  24  *      Author: Ante Bilandzic                   *
  25  *              (abilandzic@gmail.com)           *
  26  *************************************************/
  27
  28 #define AliFlowAnalysisWithCumulants_cxx
  29
  30 #include "Riostream.h"
  31 #include "TMath.h"
  32 #include "TFile.h"
  33 #include "TList.h"
  34 #include "TProfile.h"
  35 #include "TProfile2D.h"
  36 #include "TProfile3D.h"
  37 #include "TH1F.h"
  38 #include "TH1D.h"
  39 #include "TMatrixD.h"
  40 #include "TVectorD.h"
  41
  42 #include "AliFlowCommonConstants.h"
  43 #include "AliFlowCommonHist.h"
  44 #include "AliFlowCommonHistResults.h"
  45 #include "AliFlowEventSimple.h"
  46 #include "AliFlowTrackSimple.h"
  47 #include "AliFlowAnalysisWithCumulants.h"
  48 #include "AliFlowVector.h"
  49
  50 //================================================================================================================
  51
  52 ClassImp(AliFlowAnalysisWithCumulants)
  53
  54 AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants():
  55  fHistList(NULL),
  56  fHistListName(NULL),
  57  fAnalysisSettings(NULL),
  58  fCommonHists(NULL),
  59  fCommonHistsResults2nd(NULL),
  60  fCommonHistsResults4th(NULL),
  61  fCommonHistsResults6th(NULL),
  62  fCommonHistsResults8th(NULL),
  63  fnBinsPhi(0),
  64  fPhiMin(0),
  65  fPhiMax(0),
  66  fPhiBinWidth(0),
  67  fnBinsPt(0),
  68  fPtMin(0),
  69  fPtMax(0),
  70  fPtBinWidth(0),
  71  fnBinsEta(0),
  72  fEtaMin(0),
  73  fEtaMax(0),
  74  fEtaBinWidth(0),
  75  fHarmonic(2),
  76  fMultiple(1),
  77  fR0(2.2),
  78  fWeightsList(NULL),
  79  fUsePhiWeights(kFALSE),
  80  fUsePtWeights(kFALSE),
  81  fUseEtaWeights(kFALSE),
  82  fUseParticleWeights(NULL),
  83  fPhiWeights(NULL),
  84  fPtWeights(NULL),
  85  fEtaWeights(NULL),
  86  fMultiplicityWeight(NULL),
  87  fReferenceFlowList(NULL),
  88  fReferenceFlowProfiles(NULL),
  89  fReferenceFlowResults(NULL),
  90  fReferenceFlowFlags(NULL),
  91  fnBinsMult(10000),
  92  fMinMult(0.),
  93  fMaxMult(10000.),
  94  fGEBE(NULL),
  95  fReferenceFlowGenFun(NULL),
  96  fQvectorComponents(NULL),
  97  fAverageOfSquaredWeight(NULL),
  98  fAvM(0.),
  99  fnEvts(0),
 100  fReferenceFlowCumulants(NULL),
 101  fReferenceFlow(NULL),
 102  fChi(NULL),
 103  fDiffFlowList(NULL),
 104  fDiffFlowProfiles(NULL),
 105  fDiffFlowResults(NULL),
 106  fDiffFlowFlags(NULL),
 107  fTuningList(NULL),
 108  fTuningProfiles(NULL),
 109  fTuningResults(NULL),
 110  fTuningFlags(NULL),
 111  fTuneParameters(kFALSE)
 112  {
 113   // Constructor.
 114
 115   // Base list to hold all output objects:
 116   fHistList = new TList();
 117   fHistListName = new TString("cobjGFC");
 118   fHistList->SetName(fHistListName->Data());
 119   fHistList->SetOwner(kTRUE);
 120
 121   // List to hold histograms with phi, pt and eta weights:
 122   fWeightsList = new TList();
 123   fWeightsList->SetName("Weights");
 124   fWeightsList->SetOwner(kTRUE);
 125   fHistList->Add(fWeightsList);
 126
 127   // Multiplicity weight:
 128   fMultiplicityWeight = new TString("unit");
 129
 130   // Initialize all arrays:
 131   this->InitializeArrays();
 132
 133 } // end of AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants()
 134
 135 //================================================================================================================
 136
 137 AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants()
 138 {
 139  // Desctructor.
 140
 141  delete fHistList;
 142  delete fWeightsList;
 143
 144 } // end of AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants()
 145
 146 //================================================================================================================
 147
 148 void AliFlowAnalysisWithCumulants::Init()
 149 {
 150  // Initialize and book all objects.
 151
 152  // a) Cross check if the user settings make sense before starting;
 153  // b) Access all common constants;
 154  // c) Book and fill weights histograms;
 155  // d) Book and nest all lists in the base list fHistList;
 156  // e) Book and fill profile holding analysis settings;
 157  // f) Book common control and results histograms;
 158  // g) Store flags for reference flow;
 159  // h) Store flags for differential flow;
 160  // i) Book all objects needed for tuning.
 161
 162  //save old value and prevent histograms from being added to directory
 163  //to avoid name clashes in case multiple analaysis objects are used
 164  //in an analysis
 165  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 166  TH1::AddDirectory(kFALSE);
 167
 168  this->CrossCheckSettings();
 169  this->AccessConstants();
 170  this->BookAndFillWeightsHistograms();
 171  this->BookAndNestAllLists();
 172  this->BookProfileHoldingSettings();
 173  this->BookCommonHistograms();
 174  this->BookEverythingForReferenceFlow();
 175  this->BookEverythingForDiffFlow();
 176  this->StoreReferenceFlowFlags();
 177  this->StoreDiffFlowFlags();
 178  if(fTuneParameters){this->BookEverythingForTuning();}
 179
 180  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic); // to be improved (moved somewhere else?)
 181
 182  TH1::AddDirectory(oldHistAddStatus);
 183
 184 } // end of void AliFlowAnalysisWithCumulants::Init()
 185
 186 //================================================================================================================
 187
 188 void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent)
 189 {
 190  // Running over data only in this method.
 191
 192  // a) Check all pointers used in this method;
 193  // b) If tuning enabled, fill generating functions for different values of tuning parameters;
 194  // c) For default values of tuning parameters (r0 = 2.2 and cutoff at 10th order):
 195  //  c1) Fill common control histograms;
 196  //  c2) Fill generating function for reference flow;
 197  //  c3) Fill profile holding average values of various Q-vector components;
 198  //  c4) Fill generating function for differential flow.
 199
 200  this->CheckPointersUsedInMake();
 201  if(fTuneParameters) {this->FillGeneratingFunctionsForDifferentTuningParameters(anEvent);}
 202  Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of RPs (i.e. number of particles used to determine the reaction plane)
 203  if(nRP<10) {return;} // generating function formalism make sense only for nRPs >= 10 for default settings
 204  fCommonHists->FillControlHistograms(anEvent);
 205  this->FillGeneratingFunctionForReferenceFlow(anEvent);
 206  this->FillQvectorComponents(anEvent);
 207  this->FillGeneratingFunctionForDiffFlow(anEvent);
 208
 209 } // end of void AliFlowAnalysisWithCumulants::Make()
 210
 211 //================================================================================================================
 212
 213 void AliFlowAnalysisWithCumulants::Finish()
 214 {
 215  // Calculate the final results.
 216
 217  // a) Check all pointers used in this method;
 218  // b) Access settings for analysis with Generating Function Cumulants;
 219  // c) From relevant common control histogram get average multiplicity of RPs and number of events;
 220  // d) Calculate cumulants for reference flow;
 221  // e) Calculate from isotropic cumulants reference flow;
 222  // f) Calculate error for reference flow estimates;
 223  // g) Store the final results for reference flow in common hist results;
 224  // h) Print on the screen the final results for reference flow;
 225  // i) Calculate cumulants for differential flow;
 226  // j) Calculate differential flow for RPs/POIs vs pt/eta from cumulants;
 227  // k) Calculate integrated flow of RPs and POIs;
 228  // l) Print on the screen the final results for integrated flow of RPs and POIs;
 229  // m) If tuning enabled, calculate results for different tuning parameters.
 230
 231  this->CheckPointersUsedInFinish();
 232  this->AccessSettings();
 233  this->GetAvMultAndNoOfEvts();
 234  this->CalculateCumulantsForReferenceFlow();
 235  this->CalculateReferenceFlow();
 236  this->CalculateReferenceFlowError();
 237  this->FillCommonHistResultsForReferenceFlow();
 238  if(fPrintFinalResults[0]){this->PrintFinalResults("RF");}
 239  this->CalculateCumulantsForDiffFlow("RP","pt");
 240  this->CalculateCumulantsForDiffFlow("RP","eta");
 241  this->CalculateCumulantsForDiffFlow("POI","pt");
 242  this->CalculateCumulantsForDiffFlow("POI","eta");
 243  this->CalculateDifferentialFlow("RP","pt");
 244  this->CalculateDifferentialFlow("RP","eta");
 245  this->CalculateDifferentialFlow("POI","pt");
 246  this->CalculateDifferentialFlow("POI","eta");
 247  this->CalculateDifferentialFlowErrors("RP","pt");
 248  this->CalculateDifferentialFlowErrors("RP","eta");
 249  this->CalculateDifferentialFlowErrors("POI","pt");
 250  this->CalculateDifferentialFlowErrors("POI","eta");
 251  this->FillCommonHistResultsForDifferentialFlow("RP");
 252  this->FillCommonHistResultsForDifferentialFlow("POI");
 253  this->CalculateIntegratedFlow("RP");
 254  this->CalculateIntegratedFlow("POI");
 255  if(fPrintFinalResults[1]){this->PrintFinalResults("RP");}
 256  if(fPrintFinalResults[2]){this->PrintFinalResults("POI");}
 257  if(fTuneParameters){this->FinalizeTuning();}
 258
 259 } // end of void AliFlowAnalysisWithCumulants::Finish()
 260
 261 //================================================================================================================
 262
 263 void AliFlowAnalysisWithCumulants::FinalizeTuning()
 264 {
 265  // Finalize results with tuned inerpolating parameters.
 266
 267  for(Int_t r=0;r<10;r++)
 268  {
 269   if(TMath::Abs(fTuningR0[r])<1.e-10) continue; // protection against division by r0 bellow
 270   for(Int_t pq=0;pq<5;pq++)
 271   {
 272    Int_t pMax = fTuningGenFun[r][pq]->GetXaxis()->GetNbins();
 273    Int_t qMax = fTuningGenFun[r][pq]->GetYaxis()->GetNbins();
 274    // <G[p][q]>
 275    TMatrixD dAvG(pMax,qMax);
 276    dAvG.Zero();
 277    Bool_t someAvGEntryIsNegative = kFALSE;
 278    for(Int_t p=0;p<pMax;p++)
 279    {
 280     for(Int_t q=0;q<qMax;q++)
 281     {
 282      dAvG(p,q) = fTuningGenFun[r][pq]->GetBinContent(fTuningGenFun[r][pq]->GetBin(p+1,q+1));
 283      if(dAvG(p,q)<0.)
 284      {
 285       someAvGEntryIsNegative = kTRUE;
 286       cout<<endl;
 287       cout<<" WARNING: "<<Form("<G[%d][%d]> is negative !!!! GFC results are meaningless for r0 = %f, pq = %i.",p,q,fTuningR0[r],pq)<<endl;
 288       cout<<endl;
 289      }
 290     }
 291    }
 292    // C[p][q] (generating function for the cumulants)
 293    TMatrixD dC(pMax,qMax);
 294    dC.Zero();
 295    if(fAvM>0. && !someAvGEntryIsNegative)
 296    {
 297     for(Int_t p=0;p<pMax;p++)
 298     {
 299      for(Int_t q=0;q<qMax;q++)
 300      {
 301       dC(p,q) = fAvM*(pow(dAvG(p,q),(1./fAvM))-1.);
 302      }
 303     }
 304    }
 305    // Averaging the generating function for cumulants over azimuth
 306    // in order to eliminate detector effects.
 307    // <C[p][q]> (Remark: here <> stands for average over azimuth):
 308    TVectorD dAvC(pMax);
 309    dAvC.Zero();
 310    for(Int_t p=0;p<pMax;p++)
 311    {
 312     Double_t temp = 0.;
 313     for(Int_t q=0;q<qMax;q++)
 314     {
 315      temp += 1.*dC(p,q);
 316     }
 317     dAvC[p] = temp/qMax;
 318    }
 319    // Finally, the isotropic cumulants for reference flow and reference flow itself:
 320    TVectorD cumulant(pMax);
 321    cumulant.Zero();
 322    TVectorD flow(pMax);
 323    flow.Zero();
 324    if(pMax==2)
 325    {
 326     cumulant[0]=(1./(fTuningR0[r]*fTuningR0[r]))*(2.*dAvC[0]-(1./2.)*dAvC[1]);
 327     cumulant[1]=(2./pow(fTuningR0[r],4.))*((-2.)*dAvC[0]+1.*dAvC[1]);
 328     if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);}
 329     if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);}
 330    }
 331    else if(pMax==3)
 332    {
 333     cumulant[0] = (1./(fTuningR0[r]*fTuningR0[r]))*(3.*dAvC[0]-(3./2.)*dAvC[1]+(1./3.)*dAvC[2]);
 334     cumulant[1] = (2./pow(fTuningR0[r],4.))*((-5.)*dAvC[0]+4.*dAvC[1]-1.*dAvC[2]);
 335     cumulant[2] = (6./pow(fTuningR0[r],6.))*(3.*dAvC[0]-3.*dAvC[1]+1.*dAvC[2]);
 336     if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);}
 337     if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);}
 338     if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);}
 339    }
 340    else if(pMax==4)
 341    {
 342     cumulant[0] = (1./(fTuningR0[r]*fTuningR0[r]))*(4.*dAvC[0]-3.*dAvC[1]+(4./3.)*dAvC[2]-(1./4.)*dAvC[3]);
 343     cumulant[1] = (1./pow(fTuningR0[r],4.))*((-52./3.)*dAvC[0]+19.*dAvC[1]-(28./3.)*dAvC[2]+(11./6.)*dAvC[3]);
 344     cumulant[2] = (3./pow(fTuningR0[r],6.))*(18.*dAvC[0]-24.*dAvC[1]+14.*dAvC[2]-3.*dAvC[3]);
 345     cumulant[3] = (24./pow(fTuningR0[r],8.))*((-4.)*dAvC[0]+6.*dAvC[1]-4.*dAvC[2]+1.*dAvC[3]);
 346     if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);}
 347     if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);}
 348     if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);}
 349     if(cumulant[3]<=0.) {flow[3] = pow((-1./33.)*cumulant[3],1./8.);}
 350    }
 351    else if(pMax==5)
 352    {
 353     cumulant[0] = (-1./(60*fTuningR0[r]*fTuningR0[r]))*((-300.)*dAvC[0]+300.*dAvC[1]-200.*dAvC[2]+75.*dAvC[3]-12.*dAvC[4]);
 354     cumulant[1] = (-1./(6.*pow(fTuningR0[r],4.)))*(154.*dAvC[0]-214.*dAvC[1]+156.*dAvC[2]-61.*dAvC[3]+10.*dAvC[4]);
 355     cumulant[2] = (3./(2.*pow(fTuningR0[r],6.)))*(71.*dAvC[0]-118.*dAvC[1]+98.*dAvC[2]-41.*dAvC[3]+7.*dAvC[4]);
 356     cumulant[3] = (-24./pow(fTuningR0[r],8.))*(14.*dAvC[0]-26.*dAvC[1]+24.*dAvC[2]-11.*dAvC[3]+2.*dAvC[4]);
 357     cumulant[4] = (120./pow(fTuningR0[r],10.))*(5.*dAvC[0]-10.*dAvC[1]+10.*dAvC[2]-5.*dAvC[3]+1.*dAvC[4]);
 358     if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);}
 359     if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);}
 360     if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);}
 361     if(cumulant[3]<=0.) {flow[3] = pow((-1./33.)*cumulant[3],1./8.);}
 362     if(cumulant[4]>=0.) {flow[4] = pow((1./456.)*cumulant[4],1./10.);}
 363    }
 364    else if(pMax==8)
 365    {
 366     cumulant[0] = (1./(fTuningR0[r]*fTuningR0[r]))*(8.*dAvC[0]-14.*dAvC[1]+(56./3.)*dAvC[2]-(35./2.)*dAvC[3]
 367                 + (56./5.)*dAvC[4]-(14./3.)*dAvC[5]+(8./7.)*dAvC[6]-(1./8.)*dAvC[7]);
 368     cumulant[1] = (1./pow(fTuningR0[r],4.))*((-1924./35.)*dAvC[0]+(621./5.)*dAvC[1]-(8012./45.)*dAvC[2]
 369                 + (691./4.)*dAvC[3]-(564./5.)*dAvC[4]+(2143./45.)*dAvC[5]-(412./35.)*dAvC[6]+(363./280.)*dAvC[7]);
 370     cumulant[2] = (1./pow(fTuningR0[r],6.))*(349.*dAvC[0]-(18353./20.)*dAvC[1]+(7173./5.)*dAvC[2]
 371                 - 1457.*dAvC[3]+(4891./5.)*dAvC[4]-(1683./4.)*dAvC[5]+(527./5.)*dAvC[6]-(469./40.)*dAvC[7]);
 372     cumulant[3] = (1./pow(fTuningR0[r],8.))*((-10528./5.)*dAvC[0]+(30578./5.)*dAvC[1]-(51456./5.)*dAvC[2]
 373                 + 10993.*dAvC[3]-(38176./5.)*dAvC[4]+(16818./5.)*dAvC[5]-(4288./5.)*dAvC[6]+(967./10.)*dAvC[7]);
 374     cumulant[4] = (1./pow(fTuningR0[r],10.))*(11500.*dAvC[0]-35800.*dAvC[1]+63900.*dAvC[2]-71600.*dAvC[3]
 375             + 51620.*dAvC[4]-23400.*dAvC[5]+6100.*dAvC[6]-700.*dAvC[7]);
 376     cumulant[5] = (1./pow(fTuningR0[r],12.))*(-52560.*dAvC[0]+172080.*dAvC[1]-321840.*dAvC[2]+376200.*dAvC[3]
 377                 - 281520.*dAvC[4]+131760.*dAvC[5]-35280.*dAvC[6]+4140.*dAvC[7]);
 378     cumulant[6] = (1./pow(fTuningR0[r],14.))*(176400.*dAvC[0]-599760.*dAvC[1]+1164240.*dAvC[2]-1411200.*dAvC[3]
 379             + 1093680.*dAvC[4]-529200.*dAvC[5]+146160.*dAvC[6]-17640.*dAvC[7]);
 380     cumulant[7] = (1./pow(fTuningR0[r],16.))*(-322560*dAvC[0]+1128960.*dAvC[1]-2257920.*dAvC[2]+2822400.*dAvC[3]
 381                     - 2257920.*dAvC[4]+1128960.*dAvC[5]-322560.*dAvC[6]+40320.*dAvC[7]);
 382     if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);}
 383     if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);}
 384     if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);}
 385     if(cumulant[3]<=0.) {flow[3] = pow((-1./33.)*cumulant[3],1./8.);}
 386     if(cumulant[4]>=0.) {flow[4] = pow((1./456.)*cumulant[4],1./10.);}
 387     if(cumulant[5]<=0.) {flow[5] = pow((-1./9460.)*cumulant[5],1./12.);}
 388     if(cumulant[6]>=0.) {flow[6] = pow((1./274800.)*cumulant[6],1./14.);}
 389     if(cumulant[7]<=0.) {flow[7] = pow((-1./10643745.)*cumulant[7],1./16.);}
 390    }
 391    // Store cumulants and reference flow:
 392    for(Int_t co=0;co<pMax;co++) // cumulant order
 393    {
 394     fTuningCumulants[r][pq]->SetBinContent(co+1,cumulant[co]);
 395     fTuningFlow[r][pq]->SetBinContent(co+1,flow[co]);
 396    }
 397   } // end of for(Int_t pq=0;pq<5;pq++)
 398  } // end of for(Int_t r=0;r<10;r++)
 399
 400 } // end of void AliFlowAnalysisWithCumulants::FinalizeTuning()
 401
 402 //================================================================================================================
 403
 404 void AliFlowAnalysisWithCumulants::FillGeneratingFunctionsForDifferentTuningParameters(AliFlowEventSimple *anEvent)
 405 {
 406  // Fill generating function for reference flow evaluated for different tuning parameters.
 407
 408  Int_t pMax[5] = {2,3,4,5,8};
 409  Int_t qMax[5] = {5,7,9,11,17};
 410
 411  // Particle variables and weights:
 412  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 413  Double_t dPt  = 0.; // transverse momentum
 414  Double_t dEta = 0.; // pseudorapidity
 415  Double_t wPhi = 1.; // phi weight
 416  Double_t wPt  = 1.; // pt weight
 417  Double_t wEta = 1.; // eta weight
 418
 419  Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where:
 420                                           // nRP   = # of particles used to determine the reaction plane;
 421                                           // nPOI  = # of particles of interest for a detailed flow analysis;
 422                                           // rest  = # of particles which are not niether RPs nor POIs.
 423
 424  Int_t nRP = anEvent->GetEventNSelTracksRP(); // nRP = # of particles used to determine the reaction plane;
 425
 426  Double_t tuningGenFunEBE[10][5][8][17] = {{{{0.}}}};
 427  for(Int_t r=0;r<10;r++)
 428  {
 429   for(Int_t pq=0;pq<5;pq++)
 430   {
 431    for(Int_t p=0;p<pMax[pq];p++)
 432    {
 433     for(Int_t q=0;q<qMax[pq];q++)
 434     {
 435      tuningGenFunEBE[r][pq][p][q] = 1.;
 436     }
 437    }
 438   }
 439  }
 440
 441  // Looping over tracks:
 442  for(Int_t i=0;i<nPrim;i++)
 443  {
 444   AliFlowTrackSimple *aftsTrack = anEvent->GetTrack(i);
 445   if(aftsTrack && aftsTrack->InRPSelection())
 446   {
 447    // Access particle variables and weights:
 448    dPhi = aftsTrack->Phi();
 449    dPt  = aftsTrack->Pt();
 450    dEta = aftsTrack->Eta();
 451    if(fUsePhiWeights && fnBinsPhi) // determine phi weight for this particle:
 452    {
 453     wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 454    }
 455    if(fUsePtWeights && fnBinsPt) // determine pt weight for this particle:
 456    {
 457     wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 458    }
 459    if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle:
 460    {
 461     wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 462    }
 463    // Fill the generating functions:
 464    for(Int_t r=0;r<10;r++) // 10 different values for interpolating parameter r0
 465    {
 466     if(TMath::Abs(fTuningR0[r])<1.e-10) continue;
 467     for(Int_t pq=0;pq<5;pq++) // 5 different values for set (pMax,qMax)
 468     {
 469      for(Int_t p=0;p<pMax[pq];p++)
 470      {
 471       for(Int_t q=0;q<qMax[pq];q++)
 472       {
 473        tuningGenFunEBE[r][pq][p][q] *= (1.+wPhi*wPt*wEta*(2.*fTuningR0[r]*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax[pq]));
 474       } // end of for(Int_t q=0;q<qMax[pq];q++)
 475      } // end of for(Int_t p=0;p<pMax[pq];p++)
 476     } // end for(Int_t pq=0;pq<5;pq++) // 5 different values for set (pMax,qMax)
 477    } // end of for(Int_t r=0;r<10;r++) // 10 different values for interpolating parameter r0
 478   } // end of if(aftsTrack && aftsTrack->InRPSelection())
 479  } // end of for(Int_t i=0;i<nPrim;i++)
 480
 481  // Store G[p][q]:
 482  for(Int_t r=0;r<10;r++)
 483  {
 484   for(Int_t pq=0;pq<5;pq++)
 485   {
 486    for(Int_t p=0;p<pMax[pq];p++)
 487    {
 488     for(Int_t q=0;q<qMax[pq];q++)
 489     {
 490      if(fTuningGenFun[r][pq]) {fTuningGenFun[r][pq]->Fill((Double_t)p,(Double_t)q,tuningGenFunEBE[r][pq][p][q],1.);}
 491     }
 492    }
 493   }
 494  }
 495
 496 } // end of void AliFlowAnalysisWithCumulants::FillGeneratingFunctionsForDifferentTuningParameters(AliFlowEventSimple *anEvent)
 497
 498 //================================================================================================================
 499
 500 void AliFlowAnalysisWithCumulants::FillGeneratingFunctionForReferenceFlow(AliFlowEventSimple *anEvent)
 501 {
 502  // Fill generating function for reference flow for current event.
 503
 504  // Particle variables and weights:
 505  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 506  Double_t dPt  = 0.; // transverse momentum
 507  Double_t dEta = 0.; // pseudorapidity
 508  Double_t wPhi = 1.; // phi weight
 509  Double_t wPt  = 1.; // pt weight
 510  Double_t wEta = 1.; // eta weight
 511
 512  Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where:
 513                                           // nRP   = # of particles used to determine the reaction plane;
 514                                           // nPOI  = # of particles of interest for a detailed flow analysis;
 515                                           // rest  = # of particles which are not niether RPs nor POIs.
 516
 517  Int_t nRP = anEvent->GetEventNSelTracksRP(); // nRP = # of particles used to determine the reaction plane;
 518
 519  // Initializing the generating function G[p][q] for reference flow for current event:
 520  Int_t pMax = fGEBE->GetNrows();
 521  Int_t qMax = fGEBE->GetNcols();
 522  for(Int_t p=0;p<pMax;p++)
 523  {
 524   for(Int_t q=0;q<qMax;q++)
 525   {
 526    (*fGEBE)(p,q) = 1.;
 527   }
 528  }
 529
 530  // Cross-checking the number of RPs in current event:
 531  Int_t crossCheckRP = 0;
 532
 533  // Looping over tracks:
 534  for(Int_t i=0;i<nPrim;i++)
 535  {
 536   AliFlowTrackSimple *aftsTrack = anEvent->GetTrack(i);
 537   if(aftsTrack && aftsTrack->InRPSelection())
 538   {
 539    crossCheckRP++;
 540    // Access particle variables and weights:
 541    dPhi = aftsTrack->Phi();
 542    dPt  = aftsTrack->Pt();
 543    dEta = aftsTrack->Eta();
 544    if(fUsePhiWeights && fnBinsPhi) // determine phi weight for this particle:
 545    {
 546     wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 547    }
 548    if(fUsePtWeights && fnBinsPt) // determine pt weight for this particle:
 549    {
 550     wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 551    }
 552    if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle:
 553    {
 554     wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 555    }
 556    // Fill the generating function:
 557    for(Int_t p=0;p<pMax;p++)
 558    {
 559     for(Int_t q=0;q<qMax;q++)
 560     {
 561      (*fGEBE)(p,q) *= (1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax));
 562     }
 563    }
 564    // Fill the profile to calculate <<w^2>>:
 565    fAverageOfSquaredWeight->Fill(0.5,pow(wPhi*wPt*wEta,2.),1.);
 566   } // end of if(aftsTrack && aftsTrack->InRPSelection())
 567  } // end of for(Int_t i=0;i<nPrim;i++)
 568
 569  // Cross check # of RPs:
 570  if(anEvent && (crossCheckRP != anEvent->GetEventNSelTracksRP()))
 571  {
 572   cout<<endl;
 573   cout<<"WARNING (GFC): crossCheckRP != nRP in GFC::Make(). Something is wrong with RP flagging !!!!"<<endl;
 574   cout<<endl;
 575   exit(0);
 576  }
 577
 578  // Storing the value of G[p][q] in 2D profile in order to get eventually the avarage <G[p][q]>:
 579  // Determine first the event weight for G[p][q]:
 580  // (to be improved - this can be implemented much better, this shall be executed only once out of Make(), eventWeight should be a data member)
 581  Double_t eventWeight = 0.;
 582  if(!strcmp(fMultiplicityWeight->Data(),"unit"))
 583  {
 584   eventWeight = 1.;
 585  } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
 586    {
 587     eventWeight = anEvent->GetEventNSelTracksRP();
 588    }
 589  // Store G[p][q] weighted appropriately:
 590  for(Int_t p=0;p<pMax;p++)
 591  {
 592   for(Int_t q=0;q<qMax;q++)
 593   {
 594    fReferenceFlowGenFun->Fill((Double_t)p,(Double_t)q,(*fGEBE)(p,q),eventWeight);
 595   }
 596  }
 597
 598 } // end of void AliFlowAnalysisWithCumulants::FillGeneratingFunctionForReferenceFlow(AliFlowEventSimple* anEvent)
 599
 600 //================================================================================================================
 601
 602 void AliFlowAnalysisWithCumulants::FillQvectorComponents(AliFlowEventSimple* anEvent)
 603 {
 604  // Fill components of Q-vector for current event (needed for error calculation).
 605
 606  // Remark: Components are stored in profile fQvectorComponents whose binning is organized as follows:
 607  //  1st bin: Q_x
 608  //  2nd bin: Q_y
 609  //  3rd bin: (Q_x)^2
 610  //  4th bin: (Q_y)^2
 611
 612  AliFlowVector afv;
 613  afv.Set(0.,0.);
 614  afv.SetMult(0);
 615
 616  Int_t n = 2; // to be removed
 617
 618  if(anEvent)
 619  {
 620   afv = anEvent->GetQ(1*n,fWeightsList,fUsePhiWeights,fUsePtWeights,fUseEtaWeights); // get the Q-vector for this event
 621   fQvectorComponents->Fill(0.5,afv.X(),1.); // in the 1st bin fill Q_x
 622   fQvectorComponents->Fill(1.5,afv.Y(),1.); // in the 2nd bin fill Q_y
 623   fQvectorComponents->Fill(2.5,pow(afv.X(),2.),1.); // in the 3rd bin fill (Q_x)^2
 624   fQvectorComponents->Fill(3.5,pow(afv.Y(),2.),1.); // in the 4th bin fill (Q_y)^2
 625  }
 626
 627 } // end of void AliFlowAnalysisWithCumulants::FillQvectorComponents(AliFlowEventSimple* anEvent)
 628
 629 //================================================================================================================
 630
 631 void AliFlowAnalysisWithCumulants::FillGeneratingFunctionForDiffFlow(AliFlowEventSimple* anEvent)
 632 {
 633  // Fill generating function for differential flow for the current event.
 634
 635  // Remark 0: Generating function D[b][p][q] is a complex number => real and imaginary part are calculated separately
 636  //           (b denotes pt or eta bin);
 637  // Remark 1: Note that bellow G[p][q] is needed, the value of generating function for reference flow for the CURRENT event.
 638  //           This values is obtained in method FillGeneratingFunctionForReferenceFlow() as TMatrixD fGEBE;
 639  // Remark 2: Results for D[b][p][q] are stored in 3D profiles fDiffFlowGenFun[0=Re,1=Im][0=RP,1=POI][0=pt,1=eta] in order to
 640  //           automatically get <Re(D[b][p][q])> and <Im(D[b][p][q])> at the end of the day.
 641
 642  // Particle variables and weights:
 643  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
 644  Double_t dPt  = 0.; // transverse momentum
 645  Double_t dEta = 0.; // pseudorapidity
 646  Double_t wPhi = 1.; // phi weight
 647  Double_t wPt  = 1.; // pt weight
 648  Double_t wEta = 1.; // eta weight
 649
 650  // pMax and qMax:
 651  Int_t pMax = fGEBE->GetNrows();
 652  Int_t qMax = fGEBE->GetNcols();
 653
 654  Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where:
 655                                           // nRP   = # of particles used to determine the reaction plane;
 656                                           // nPOI  = # of particles of interest for a detailed flow analysis;
 657                                           // rest  = # of particles which are not niether RPs nor POIs.
 658
 659  Int_t nRP = anEvent->GetEventNSelTracksRP(); // nRP = # of particles used to determine the reaction plane
 660
 661  // Start the second loop over event in order to evaluate the generating function D[b][p][q] for differential flow:
 662  for(Int_t i=0;i<nPrim;i++)
 663  {
 664   AliFlowTrackSimple *aftsTrack = anEvent->GetTrack(i);
 665   if(aftsTrack)
 666   {
 667    if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())) continue;
 668    // Differential flow of POIs:
 669    if(aftsTrack->InPOISelection())
 670    {
 671     // Get azimuthal angle, momentum and pseudorapidity of a particle:
 672     dPhi = aftsTrack->Phi();
 673     dPt  = aftsTrack->Pt();
 674     dEta = aftsTrack->Eta();
 675     Double_t ptEta[2] = {dPt,dEta};
 676
 677     // Count number of POIs in pt/eta bin:
 678     for(Int_t pe=0;pe<2;pe++)
 679     {
 680      fNoOfParticlesInBin[1][pe]->Fill(ptEta[pe],ptEta[pe],1.);
 681     }
 682
 683     if(!(aftsTrack->InRPSelection())) // particle was flagged only as POI
 684     {
 685      // Fill generating function:
 686      for(Int_t p=0;p<pMax;p++)
 687      {
 688       for(Int_t q=0;q<qMax;q++)
 689       {
 690        for(Int_t ri=0;ri<2;ri++)
 691        {
 692         for(Int_t pe=0;pe<2;pe++)
 693         {
 694          if(ri==0) // Real part (to be improved - this can be implemented better)
 695          {
 696           fDiffFlowGenFun[ri][1][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow
 697                                        (*fGEBE)(p,q)*cos(fMultiple*fHarmonic*dPhi),1.);
 698          }
 699          else if(ri==1) // Imaginary part (to be improved - this can be implemented better)
 700          {
 701           fDiffFlowGenFun[ri][1][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow
 702                                        (*fGEBE)(p,q)*sin(fMultiple*fHarmonic*dPhi),1.);
 703          }
 704         } // end of for(Int_t pe=0;pe<2;pe++)
 705        } // end of for(Int_t ri=0;ri<2;ri++)
 706       } // end of for(Int_t q=0;q<qMax;q++)
 707      } // end of for(Int_t p=0;p<pMax;p++)
 708     } // end of if(!(aftsTrack->InRPSelection())) // particle was flagged only as POI
 709     else if(aftsTrack->InRPSelection()) // particle was flagged both as RP and POI
 710     {
 711      // If particle weights were used, get them:
 712      if(fUsePhiWeights && fnBinsPhi) // determine phi weight for this particle:
 713      {
 714       wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 715      }
 716      if(fUsePtWeights && fnBinsPt) // determine pt weight for this particle:
 717      {
 718       wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 719      }
 720      if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle:
 721      {
 722       wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 723      }
 724      // Fill generating function:
 725      for(Int_t p=0;p<pMax;p++)
 726      {
 727       for(Int_t q=0;q<qMax;q++)
 728       {
 729        for(Int_t ri=0;ri<2;ri++)
 730        {
 731         for(Int_t pe=0;pe<2;pe++)
 732         {
 733          if(ri==0) // Real part (to be improved - this can be implemented better)
 734          {
 735           fDiffFlowGenFun[ri][1][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow
 736                                       (*fGEBE)(p,q)*cos(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.);
 737          }
 738          else if(ri==1) // Imaginary part (to be improved - this can be implemented better)
 739          {
 740           fDiffFlowGenFun[ri][1][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow
 741                                       (*fGEBE)(p,q)*sin(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.);
 742          }
 743         } // end of for(Int_t pe=0;pe<2;pe++)
 744        } // end of for(Int_t ri=0;ri<2;ri++)
 745       } // end of for(Int_t q=0;q<qMax;q++)
 746      } // end of for(Int_t p=0;p<pMax;p++)
 747     } // end of else if (aftsTrack->InRPSelection()) // particle was flagged both as RP and POI
 748    } // end of if(aftsTrack->InPOISelection())
 749    // Differential flow of RPs:
 750    if(aftsTrack->InRPSelection())
 751    {
 752     // Get azimuthal angle, momentum and pseudorapidity of a particle:
 753     dPhi = aftsTrack->Phi();
 754     dPt  = aftsTrack->Pt();
 755     dEta = aftsTrack->Eta();
 756     Double_t ptEta[2] = {dPt,dEta};
 757
 758     // Count number of RPs in pt/eta bin:
 759     for(Int_t pe=0;pe<2;pe++)
 760     {
 761      fNoOfParticlesInBin[0][pe]->Fill(ptEta[pe],ptEta[pe],1.);
 762     }
 763
 764     // If particle weights were used, get them:
 765     if(fUsePhiWeights && fnBinsPhi) // determine phi weight for this particle:
 766     {
 767      wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
 768     }
 769     if(fUsePtWeights && fnBinsPt) // determine pt weight for this particle:
 770     {
 771      wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth)));
 772     }
 773     if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle:
 774     {
 775      wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth)));
 776     }
 777     // Fill generating function:
 778     for(Int_t p=0;p<pMax;p++)
 779     {
 780      for(Int_t q=0;q<qMax;q++)
 781      {
 782       for(Int_t ri=0;ri<2;ri++)
 783       {
 784        for(Int_t pe=0;pe<2;pe++)
 785        {
 786         if(ri==0) // Real part (to be improved - this can be implemented better)
 787         {
 788          fDiffFlowGenFun[ri][0][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow
 789                                      (*fGEBE)(p,q)*cos(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.);
 790         }
 791         else if(ri==1) // Imaginary part (to be improved - this can be implemented better)
 792         {
 793          fDiffFlowGenFun[ri][0][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow
 794                                      (*fGEBE)(p,q)*sin(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.);
 795         }
 796        } // end of for(Int_t pe=0;pe<2;pe++)
 797       } // end of for(Int_t ri=0;ri<2;ri++)
 798      } // end of for(Int_t q=0;q<qMax;q++)
 799     } // end of for(Int_t p=0;p<pMax;p++)
 800    } // end of if(aftsTrack->InRPSelection())
 801   } // end of if(aftsTrack)
 802  } // end of for(Int_t i=0;i<nPrim;i++)
 803
 804 } // end of void AliFlowAnalysisWithCumulants::FillGeneratingFunctionForDiffFlow(AliFlowEventSimple* anEvent)
 805
 806 //================================================================================================================
 807
 808 void AliFlowAnalysisWithCumulants::GetOutputHistograms(TList *outputListHistos)
 809 {
 810  // Get pointers to all objects saved in the output file.
 811
 812  if(outputListHistos)
 813  {
 814   this->SetHistList(outputListHistos);
 815   if(!fHistList)
 816   {
 817    cout<<endl;
 818    cout<<" WARNING (GFC): fHistList is NULL in AFAWGFC::GOH() !!!!"<<endl;
 819    cout<<endl;
 820    exit(0);
 821   }
 822   this->GetPointersForBaseHistograms();
 823   this->GetPointersForCommonControlHistograms();
 824   this->GetPointersForCommonResultsHistograms();
 825   this->GetPointersForParticleWeightsHistograms();
 826   this->GetPointersForReferenceFlowObjects();
 827   this->GetPointersForDiffFlowObjects();
 828  } else
 829    {
 830     cout<<endl;
 831     cout<<" WARNING (GFC): outputListHistos is NULL in AFAWGFC::GOH() !!!!"<<endl;
 832     cout<<endl;
 833     exit(0);
 834    }
 835
 836 } // end of void AliFlowAnalysisWithCumulants::GetOutputHistograms(TList *outputListHistos)
 837
 838 //================================================================================================================
 839
 840 void AliFlowAnalysisWithCumulants::GetPointersForBaseHistograms()
 841 {
 842  // Get pointers to base histograms.
 843
 844  TString analysisSettingsName = "fAnalysisSettings";
 845  TProfile *analysisSettings = dynamic_cast<TProfile*>(fHistList->FindObject(analysisSettingsName.Data()));
 846  if(analysisSettings)
 847  {
 848   this->SetAnalysisSettings(analysisSettings);
 849  } else
 850    {
 851     cout<<endl;
 852     cout<<" WARNING (GFC): analysisSettings is NULL in AFAWGFC::GPFBH() !!!!"<<endl;
 853     cout<<endl;
 854     exit(0);
 855    }
 856
 857 } // end of void AliFlowAnalysisWithCumulants::GetPointersForBaseHistograms()
 858
 859 //================================================================================================================
 860
 861 void AliFlowAnalysisWithCumulants::GetPointersForCommonControlHistograms()
 862 {
 863  // Get pointers for common control histograms.
 864
 865  TString commonHistsName = "AliFlowCommonHistGFC";
 866  AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(fHistList->FindObject(commonHistsName.Data()));
 867  if(commonHist)
 868  {
 869   this->SetCommonHists(commonHist);
 870  } else
 871    {
 872     cout<<endl;
 873     cout<<" WARNING (GFC): commonHist is NULL in AFAWGFC::GPFCH() !!!!"<<endl;
 874     cout<<endl;
 875     exit(0);
 876    }
 877
 878 } // end of void AliFlowAnalysisWithCumulants::GetPointersForCommonControlHistograms()
 879
 880 //================================================================================================================
 881
 882 void AliFlowAnalysisWithCumulants::GetPointersForCommonResultsHistograms()
 883 {
 884  // Get pointers for common results histograms.
 885
 886  TString commonHistResults2ndOrderName = "AliFlowCommonHistResults2ndOrderGFC";
 887  AliFlowCommonHistResults *commonHistRes2nd = dynamic_cast<AliFlowCommonHistResults*>
 888                                               (fHistList->FindObject(commonHistResults2ndOrderName.Data()));
 889  if(commonHistRes2nd)
 890  {
 891   this->SetCommonHistsResults2nd(commonHistRes2nd);
 892  } else
 893    {
 894     cout<<endl;
 895     cout<<" WARNING (GFC): commonHistRes2nd is NULL in AFAWGFC::GPFCRH() !!!!"<<endl;
 896     cout<<endl;
 897     exit(0);
 898    }
 899  TString commonHistResults4thOrderName = "AliFlowCommonHistResults4thOrderGFC";
 900  AliFlowCommonHistResults *commonHistRes4th = dynamic_cast<AliFlowCommonHistResults*>
 901                                               (fHistList->FindObject(commonHistResults4thOrderName.Data()));
 902  if(commonHistRes4th)
 903  {
 904   this->SetCommonHistsResults4th(commonHistRes4th);
 905  } else
 906    {
 907     cout<<endl;
 908     cout<<" WARNING (GFC): commonHistRes4th is NULL in AFAWGFC::GPFCRH() !!!!"<<endl;
 909     cout<<endl;
 910     exit(0);
 911    }
 912  TString commonHistResults6thOrderName = "AliFlowCommonHistResults6thOrderGFC";
 913  AliFlowCommonHistResults *commonHistRes6th = dynamic_cast<AliFlowCommonHistResults*>
 914                                               (fHistList->FindObject(commonHistResults6thOrderName.Data()));
 915  if(commonHistRes6th)
 916  {
 917   this->SetCommonHistsResults6th(commonHistRes6th);
 918  } else
 919    {
 920     cout<<endl;
 921     cout<<" WARNING (GFC): commonHistRes6th is NULL in AFAWGFC::GPFCRH() !!!!"<<endl;
 922     cout<<endl;
 923     exit(0);
 924    }
 925  TString commonHistResults8thOrderName = "AliFlowCommonHistResults8thOrderGFC";
 926  AliFlowCommonHistResults *commonHistRes8th = dynamic_cast<AliFlowCommonHistResults*>
 927                                               (fHistList->FindObject(commonHistResults8thOrderName.Data()));
 928  if(commonHistRes8th)
 929  {
 930   this->SetCommonHistsResults8th(commonHistRes8th);
 931  } else
 932    {
 933     cout<<endl;
 934     cout<<" WARNING (GFC): commonHistRes8th is NULL in AFAWGFC::GPFCRH() !!!!"<<endl;
 935     cout<<endl;
 936     exit(0);
 937    }
 938
 939 } // end of void AliFlowAnalysisWithCumulants::GetPointersForCommonResultsHistograms()
 940
 941 //================================================================================================================
 942
 943 void AliFlowAnalysisWithCumulants::GetPointersForParticleWeightsHistograms()
 944 {
 945  // Get pointers for histograms holding particle weights.
 946
 947  TList *weightsList = dynamic_cast<TList*>(fHistList->FindObject("Weights"));
 948  if(weightsList) this->SetWeightsList(weightsList);
 949  TString fUseParticleWeightsName = "fUseParticleWeights";
 950  TProfile *useParticleWeights = dynamic_cast<TProfile*>(weightsList->FindObject(fUseParticleWeightsName.Data()));
 951  if(useParticleWeights)
 952  {
 953   this->SetUseParticleWeights(useParticleWeights);
 954   fUsePhiWeights = (Int_t)fUseParticleWeights->GetBinContent(1);
 955   fUsePtWeights = (Int_t)fUseParticleWeights->GetBinContent(2);
 956   fUseEtaWeights = (Int_t)fUseParticleWeights->GetBinContent(3);
 957  }
 958
 959 } // end of void AliFlowAnalysisWithCumulants::GetPointersForParticleWeightsHistograms();
 960
 961 //================================================================================================================
 962
 963 void AliFlowAnalysisWithCumulants::GetPointersForReferenceFlowObjects()
 964 {
 965  // Get pointers for all objects relevant for calculation of reference flow.
 966
 967  // a) Get pointers to all lists relevant for reference flow;
 968  // b) Get pointer to profile holding flags;
 969  // c) Get pointers to all objects in the list fReferenceFlowProfiles;
 970  // d) Get pointers to all objects in the list fReferenceFlowResults.
 971
 972  // a) Get pointers to all lists relevant for reference flow:
 973  TList *referenceFlowList = dynamic_cast<TList*>(fHistList->FindObject("Reference Flow"));
 974  if(!referenceFlowList)
 975  {
 976   cout<<endl;
 977   cout<<"WARNING (GFC): referenceFlowList is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
 978   cout<<endl;
 979   exit(0);
 980  }
 981  TList *referenceFlowProfiles = dynamic_cast<TList*>(referenceFlowList->FindObject("Profiles"));
 982  if(!referenceFlowProfiles)
 983  {
 984   cout<<endl;
 985   cout<<"WARNING (GFC): referenceFlowProfiles is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
 986   cout<<endl;
 987   exit(0);
 988  }
 989  TList *referenceFlowResults = dynamic_cast<TList*>(referenceFlowList->FindObject("Results"));
 990  if(!referenceFlowResults)
 991  {
 992   cout<<endl;
 993   cout<<"WARNING (GFC): referenceFlowResults is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
 994   cout<<endl;
 995   exit(0);
 996  }
 997
 998  // b) Get pointer to profile holding flags:
 999  TString referenceFlowFlagsName = "fReferenceFlowFlags";
1000  TProfile *referenceFlowFlags = dynamic_cast<TProfile*>(referenceFlowList->FindObject(referenceFlowFlagsName.Data()));
1001  if(referenceFlowFlags)
1002  {
1003   this->SetReferenceFlowFlags(referenceFlowFlags);
1004  } else
1005    {
1006     cout<<endl;
1007     cout<<"WARNING (GFC): referenceFlowFlags is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
1008     cout<<endl;
1009     exit(0);
1010    }
1011
1012  // c) Get pointers to all objects in the list fReferenceFlowProfiles:
1013  TString referenceFlowGenFunName = "fReferenceFlowGenFun";
1014  TProfile2D *referenceFlowGenFun = dynamic_cast<TProfile2D*>(referenceFlowProfiles->FindObject(referenceFlowGenFunName.Data()));
1015  if(referenceFlowGenFun)
1016  {
1017   this->SetReferenceFlowGenFun(referenceFlowGenFun);
1018  } else
1019    {
1020     cout<<endl;
1021     cout<<"WARNING (GFC): referenceFlowGenFun is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
1022     cout<<endl;
1023     exit(0);
1024    }
1025  TString qvectorComponentsName = "fQvectorComponents";
1026  TProfile *qvectorComponents = dynamic_cast<TProfile*>(referenceFlowProfiles->FindObject(qvectorComponentsName.Data()));
1027  if(qvectorComponents)
1028  {
1029   this->SetQvectorComponents(qvectorComponents);
1030  } else
1031    {
1032     cout<<endl;
1033     cout<<"WARNING (GFC): qvectorComponents is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
1034     cout<<endl;
1035     exit(0);
1036    }
1037  // <<w^2>>, where w = wPhi*wPt*wEta:
1038  TString averageOfSquaredWeightName = "fAverageOfSquaredWeight";
1039  TProfile *averageOfSquaredWeight = dynamic_cast<TProfile*>(referenceFlowProfiles->FindObject(averageOfSquaredWeightName.Data()));
1040  if(averageOfSquaredWeight)
1041  {
1042   this->SetAverageOfSquaredWeight(averageOfSquaredWeight);
1043  } else
1044    {
1045     cout<<endl;
1046     cout<<"WARNING (GFC): averageOfSquaredWeight is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
1047     cout<<endl;
1048     exit(0);
1049    }
1050
1051  // d) Get pointers to all objects in the list fReferenceFlowResults:
1052  // Final results for isotropic cumulants for reference flow:
1053  TString referenceFlowCumulantsName = "fReferenceFlowCumulants";
1054  TH1D *referenceFlowCumulants = dynamic_cast<TH1D*>(referenceFlowResults->FindObject(referenceFlowCumulantsName.Data()));
1055  if(referenceFlowCumulants)
1056  {
1057   this->SetReferenceFlowCumulants(referenceFlowCumulants);
1058  } else
1059    {
1060     cout<<endl;
1061     cout<<"WARNING (GFC): referenceFlowCumulants is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
1062     cout<<endl;
1063     exit(0);
1064    }
1065  // Final results for reference flow:
1066  TString referenceFlowName = "fReferenceFlow";
1067  TH1D *referenceFlow = dynamic_cast<TH1D*>(referenceFlowResults->FindObject(referenceFlowName.Data()));
1068  if(referenceFlow)
1069  {
1070   this->SetReferenceFlow(referenceFlow);
1071  } else
1072    {
1073     cout<<endl;
1074     cout<<"WARNING (GFC): referenceFlow is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
1075     cout<<endl;
1076     exit(0);
1077    }
1078  // Final results for resolution:
1079  TString chiName = "fChi";
1080  TH1D *chi = dynamic_cast<TH1D*>(referenceFlowResults->FindObject(chiName.Data()));
1081  if(chi)
1082  {
1083   this->SetChi(chi);
1084  } else
1085    {
1086     cout<<endl;
1087     cout<<"WARNING (GFC): chi is NULL in AFAWGFC::GPFRFO() !!!!"<<endl;
1088     cout<<endl;
1089     exit(0);
1090    }
1091
1092 } // end of void AliFlowAnalysisWithCumulants::GetPointersForReferenceFlowObjects()
1093
1094 //================================================================================================================
1095
1096 void AliFlowAnalysisWithCumulants::GetPointersForDiffFlowObjects()
1097 {
1098  // Get pointers to all objects relevant for differential flow.
1099
1100  //  a) Define flags locally (to be improved: should I promote flags to data members?);
1101  //  b) Get pointer to base list for differential flow fDiffFlowList and nested lists fDiffFlowListProfiles and fDiffFlowListResults;
1102  //  c) Get pointer to profile fDiffFlowFlags holding all flags for differential flow;
1103  //  d) Get pointers to all profiles in the list fDiffFlowProfiles;
1104  //  e) Get pointers to all profiles in the list fDiffFlowResults.
1105
1106  // a) Define flags locally (to be improved: should I promote flags to data members?):
1107  TString reIm[2] = {"Re","Im"};
1108  TString rpPoi[2] = {"RP","POI"};
1109  TString ptEta[2] = {"p_{t}","#eta"};
1110  TString order[4] = {"2nd order","4th order","6th order","8th order"};
1111  //TString differentialFlowIndex[4] = {"v'{2}","v'{4}","v'{6}","v'{8}"};
1112
1113  // b) Get pointer to base list for differential flow fDiffFlowList and nested lists fDiffFlowListProfiles and fDiffFlowListResults:
1114  TList *diffFlowList = dynamic_cast<TList*>(fHistList->FindObject("Differential Flow")); // to be improved (hardwired name)
1115  if(!diffFlowList)
1116  {
1117   cout<<endl;
1118   cout<<"WARNING: diffFlowList is NULL in AFAWQC::GPFDFO() !!!!"<<endl;
1119   cout<<endl;
1120   exit(0);
1121  }
1122  TList *diffFlowProfiles = dynamic_cast<TList*>(diffFlowList->FindObject("Profiles")); // to be improved (hardwired name)
1123  if(!diffFlowProfiles)
1124  {
1125   cout<<endl;
1126   cout<<"WARNING: diffFlowProfiles is NULL in AFAWQC::GPFDFO() !!!!"<<endl;
1127   cout<<endl;
1128   exit(0);
1129  }
1130  TList *diffFlowResults = dynamic_cast<TList*>(diffFlowList->FindObject("Results")); // to be improved (hardwired name)
1131  if(!diffFlowResults)
1132  {
1133   cout<<endl;
1134   cout<<"WARNING: diffFlowResults is NULL in AFAWQC::GPFDFO() !!!!"<<endl;
1135   cout<<endl;
1136   exit(0);
1137  }
1138
1139  // c) Get pointer to profile holding flags:
1140  TString diffFlowFlagsName = "fDiffFlowFlags";
1141  TProfile *diffFlowFlags = dynamic_cast<TProfile*>(diffFlowList->FindObject(diffFlowFlagsName.Data()));
1142  if(diffFlowFlags)
1143  {
1144   this->SetDiffFlowFlags(diffFlowFlags);
1145  } else
1146    {
1147     cout<<endl;
1148     cout<<"WARNING (GFC): diffFlowFlags is NULL in AFAWGFC::GPFDFO() !!!!"<<endl;
1149     cout<<endl;
1150     exit(0);
1151    }
1152
1153  // d) Get pointers to all profiles in the list fDiffFlowListProfiles:
1154  // Generating functions for differential flow:
1155  TProfile3D *diffFlowGenFun[2][2][2] = {{{NULL}}};
1156  for(Int_t ri=0;ri<2;ri++)
1157  {
1158   for(Int_t rp=0;rp<2;rp++)
1159   {
1160    for(Int_t pe=0;pe<2;pe++)
1161    {
1162     diffFlowGenFun[ri][rp][pe] = dynamic_cast<TProfile3D*> // to be improved - harwired name fDiffFlowGenFun in the line bellow
1163                                  (diffFlowProfiles->FindObject(Form("fDiffFlowGenFun (%s, %s, %s)",reIm[ri].Data(),rpPoi[rp].Data(),ptEta[pe].Data())));
1164     if(diffFlowGenFun[ri][rp][pe])
1165     {
1166      this->SetDiffFlowGenFun(diffFlowGenFun[ri][rp][pe],ri,rp,pe);
1167     } else
1168       {
1169        cout<<endl;
1170        cout<<"WARNING (GFC): "<<Form("diffFlowGenFun[%d][%d][%d]",ri,rp,pe)<<" is NULL in AFAWGFC::GPFDFO() !!!!"<<endl;
1171        cout<<endl;
1172        exit(0);
1173       }
1174    }
1175   }
1176  }
1177  // Number of particles in pt/eta bin for RPs/POIs:
1178  TProfile *noOfParticlesInBin[2][2] = {{NULL}};
1179  for(Int_t rp=0;rp<2;rp++)
1180  {
1181   for(Int_t pe=0;pe<2;pe++)
1182   {
1183    noOfParticlesInBin[rp][pe] = dynamic_cast<TProfile*> // to be improved - harwired name fNoOfParticlesInBin in the line bellow
1184                                  (diffFlowProfiles->FindObject(Form("fNoOfParticlesInBin (%s, %s)",rpPoi[rp].Data(),ptEta[pe].Data())));
1185    if(noOfParticlesInBin[rp][pe])
1186    {
1187     this->SetNoOfParticlesInBin(noOfParticlesInBin[rp][pe],rp,pe);
1188    } else
1189      {
1190       cout<<endl;
1191       cout<<"WARNING (GFC): "<<Form("noOfParticlesInBin[%d][%d]",rp,pe)<<" is NULL in AFAWGFC::GPFDFO() !!!!"<<endl;
1192       cout<<endl;
1193       exit(0);
1194      }
1195   }
1196  }
1197  // Differential cumulants per pt/eta bin for RPs/POIs:
1198  TH1D *diffFlowCumulants[2][2][4] = {{{NULL}}};
1199  for(Int_t rp=0;rp<2;rp++)
1200  {
1201   for(Int_t pe=0;pe<2;pe++)
1202   {
1203    for(Int_t co=0;co<4;co++)
1204    {
1205     diffFlowCumulants[rp][pe][co] = dynamic_cast<TH1D*> // to be improved - hardwired name fDiffFlowCumulants in the line bellow
1206                                     (diffFlowResults->FindObject(Form("fDiffFlowCumulants (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data())));
1207     if(diffFlowCumulants[rp][pe][co])
1208     {
1209      this->SetDiffFlowCumulants(diffFlowCumulants[rp][pe][co],rp,pe,co);
1210     } else
1211       {
1212        cout<<endl;
1213        cout<<"WARNING (GFC): "<<Form("diffFlowCumulants[%d][%d][%d]",rp,pe,co)<<" is NULL in AFAWGFC::GPFDFO() !!!!"<<endl;
1214        cout<<endl;
1215        exit(0);
1216       }
1217    }
1218   }
1219  }
1220  // Differential flow per pt/eta bin for RPs/POIs:
1221  TH1D *diffFlow[2][2][4] = {{{NULL}}};
1222  for(Int_t rp=0;rp<2;rp++)
1223  {
1224   for(Int_t pe=0;pe<2;pe++)
1225   {
1226    for(Int_t co=0;co<4;co++)
1227    {
1228     diffFlow[rp][pe][co] = dynamic_cast<TH1D*> // to be improved - hardwired name fDiffFlow in the line bellow
1229                            (diffFlowResults->FindObject(Form("fDiffFlow (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data())));
1230     if(diffFlow[rp][pe][co])
1231     {
1232      this->SetDiffFlow(diffFlow[rp][pe][co],rp,pe,co);
1233     } else
1234       {
1235        cout<<endl;
1236        cout<<"WARNING (GFC): "<<Form("diffFlow[%d][%d][%d]",rp,pe,co)<<" is NULL in AFAWGFC::GPFDFO() !!!!"<<endl;
1237        cout<<endl;
1238        exit(0);
1239       }
1240    }
1241   }
1242  }
1243
1244 } // end of void AliFlowAnalysisWithCumulants::GetPointersForDiffFlowObjects()
1245
1246 //================================================================================================================
1247
1248 void AliFlowAnalysisWithCumulants::CalculateIntegratedFlow(TString rpPoi)
1249 {
1250  // Calculate final results for integrated flow of RPs and POIs.
1251  // (to be improved - this method can be implemented much better)
1252
1253  Int_t rp = -1;
1254
1255  if(rpPoi == "RP")
1256  {
1257   rp = 0;
1258  } else if(rpPoi == "POI")
1259    {
1260     rp = 1;
1261    }
1262
1263  // pt yield:
1264  TH1F *yieldPt = NULL;
1265
1266  if(rpPoi == "POI")
1267  {
1268   yieldPt = (TH1F*)(fCommonHists->GetHistPtPOI())->Clone();
1269  } else if(rpPoi == "RP")
1270    {
1271     yieldPt = (TH1F*)(fCommonHists->GetHistPtRP())->Clone();
1272    }
1273
1274  TH1D *flow2ndPt = (TH1D*)fDiffFlow[rp][0][0]->Clone();
1275  TH1D *flow4thPt = (TH1D*)fDiffFlow[rp][0][1]->Clone();
1276  TH1D *flow6thPt = (TH1D*)fDiffFlow[rp][0][2]->Clone();
1277  TH1D *flow8thPt = (TH1D*)fDiffFlow[rp][0][3]->Clone();
1278  Double_t dvn2nd = 0., dvn4th = 0., dvn6th = 0., dvn8th = 0.; // differential flow
1279  Double_t dErrvn2nd = 0., dErrvn4th = 0., dErrvn6th = 0., dErrvn8th = 0.; // error on differential flow
1280  Double_t dVn2nd = 0., dVn4th = 0., dVn6th = 0., dVn8th = 0.; // integrated flow
1281  Double_t dErrVn2nd = 0., dErrVn4th = 0., dErrVn6th = 0., dErrVn8th = 0.; // error on integrated flow
1282  Double_t dYield = 0.; // pt yield
1283  Double_t dSum2nd = 0., dSum4th = 0., dSum6th = 0., dSum8th = 0.; // needed for normalizing integrated flow
1284  fnBinsPt = flow2ndPt->GetXaxis()->GetNbins();
1285  // looping over pt bins:
1286  for(Int_t p=1;p<=fnBinsPt;p++)
1287  {
1288   dvn2nd = flow2ndPt->GetBinContent(p);
1289   dvn4th = flow4thPt->GetBinContent(p);
1290   dvn6th = flow6thPt->GetBinContent(p);
1291   dvn8th = flow8thPt->GetBinContent(p);
1292
1293   dErrvn2nd = flow2ndPt->GetBinError(p);
1294   dErrvn4th = flow4thPt->GetBinError(p);
1295   dErrvn6th = flow6thPt->GetBinError(p);
1296   dErrvn8th = flow8thPt->GetBinError(p);
1297
1298   dYield = yieldPt->GetBinContent(p);
1299
1300   dVn2nd += dvn2nd*dYield;
1301   dVn4th += dvn4th*dYield;
1302   dVn6th += dvn6th*dYield;
1303   dVn8th += dvn8th*dYield;
1304
1305   dSum2nd += dYield;
1306   dSum4th += dYield;
1307   dSum6th += dYield;
1308   dSum8th += dYield;
1309
1310   dErrVn2nd += dYield*dYield*dErrvn2nd*dErrvn2nd;
1311   dErrVn4th += dYield*dYield*dErrvn4th*dErrvn4th;
1312   dErrVn6th += dYield*dYield*dErrvn6th*dErrvn6th;
1313   dErrVn8th += dYield*dYield*dErrvn8th*dErrvn8th;
1314  } // end of for(Int_t p=1;p<=fnBinsPt;p++)
1315
1316  // normalizing the results for integrated flow:
1317  if(dSum2nd)
1318  {
1319   dVn2nd /= dSum2nd;
1320   dErrVn2nd /= (dSum2nd*dSum2nd);
1321   dErrVn2nd = TMath::Sqrt(dErrVn2nd);
1322  }
1323  if(dSum4th)
1324  {
1325   dVn4th /= dSum4th;
1326   dErrVn4th /= (dSum4th*dSum4th);
1327   dErrVn4th = TMath::Sqrt(dErrVn4th);
1328  }
1329  if(dSum6th)
1330  {
1331   dVn6th /= dSum6th;
1332   dErrVn6th /= (dSum6th*dSum6th);
1333   dErrVn6th = TMath::Sqrt(dErrVn6th);
1334  }
1335  if(dSum8th)
1336  {
1337   dVn8th /= dSum8th;
1338   dErrVn8th /= (dSum8th*dSum8th);
1339   dErrVn8th = TMath::Sqrt(dErrVn8th);
1340  }
1341
1342  // storing the results for integrated flow in common hist results:
1343  if(rpPoi == "POI")
1344  {
1345   fCommonHistsResults2nd->FillIntegratedFlowPOI(dVn2nd,dErrVn2nd);
1346   fCommonHistsResults4th->FillIntegratedFlowPOI(dVn4th,dErrVn4th);
1347   fCommonHistsResults6th->FillIntegratedFlowPOI(dVn6th,dErrVn6th);
1348   fCommonHistsResults8th->FillIntegratedFlowPOI(dVn8th,dErrVn8th);
1349  }
1350  else if(rpPoi == "RP")
1351  {
1352   fCommonHistsResults2nd->FillIntegratedFlowRP(dVn2nd,dErrVn2nd);
1353   fCommonHistsResults4th->FillIntegratedFlowRP(dVn4th,dErrVn4th);
1354   fCommonHistsResults6th->FillIntegratedFlowRP(dVn6th,dErrVn6th);
1355   fCommonHistsResults8th->FillIntegratedFlowRP(dVn8th,dErrVn8th);
1356  }
1357
1358  delete flow2ndPt;
1359  delete flow4thPt;
1360  delete flow6thPt;
1361  delete flow8thPt;
1362  delete yieldPt;
1363
1364 } // end of void AliFlowAnalysisWithCumulants::CalculateIntegratedFlow(TString rpPoi)
1365
1366 //================================================================================================================
1367
1368 void AliFlowAnalysisWithCumulants::FillCommonHistResultsForDifferentialFlow(TString rpPoi)
1369 {
1370  // Fill common result histograms for differential flow.
1371  // (to be improved - this method can be implemented much better)
1372
1373  Int_t rp = -1;
1374
1375  if(rpPoi == "RP")
1376  {
1377   rp = 0;
1378  } else if(rpPoi == "POI")
1379    {
1380     rp = 1;
1381    }
1382
1383  // pt:
1384  for(Int_t p=1;p<=fnBinsPt;p++)
1385  {
1386   // Result:
1387   Double_t v2 = fDiffFlow[rp][0][0]->GetBinContent(p);
1388   Double_t v4 = fDiffFlow[rp][0][1]->GetBinContent(p);
1389   Double_t v6 = fDiffFlow[rp][0][2]->GetBinContent(p);
1390   Double_t v8 = fDiffFlow[rp][0][3]->GetBinContent(p);
1391   // Error:
1392   Double_t v2Error = fDiffFlow[rp][0][0]->GetBinError(p);
1393   Double_t v4Error = fDiffFlow[rp][0][1]->GetBinError(p);
1394   Double_t v6Error = fDiffFlow[rp][0][2]->GetBinError(p);
1395   Double_t v8Error = fDiffFlow[rp][0][3]->GetBinError(p);
1396   // Fill common hist results:
1397   if(rpPoi == "RP")
1398   {
1399    fCommonHistsResults2nd->FillDifferentialFlowPtRP(p,v2,v2Error);
1400    fCommonHistsResults4th->FillDifferentialFlowPtRP(p,v4,v4Error);
1401    fCommonHistsResults6th->FillDifferentialFlowPtRP(p,v6,v6Error);
1402    fCommonHistsResults8th->FillDifferentialFlowPtRP(p,v8,v8Error);
1403   } else if(rpPoi == "POI")
1404     {
1405      fCommonHistsResults2nd->FillDifferentialFlowPtPOI(p,v2,v2Error);
1406      fCommonHistsResults4th->FillDifferentialFlowPtPOI(p,v4,v4Error);
1407      fCommonHistsResults6th->FillDifferentialFlowPtPOI(p,v6,v6Error);
1408      fCommonHistsResults8th->FillDifferentialFlowPtPOI(p,v8,v8Error);
1409     }
1410  } // end of for(Int_t p=1;p<=fnBinsPt;p++)
1411
1412  // eta:
1413  for(Int_t e=1;e<=fnBinsEta;e++)
1414  {
1415   // Results:
1416   Double_t v2 = fDiffFlow[rp][1][0]->GetBinContent(e);
1417   Double_t v4 = fDiffFlow[rp][1][1]->GetBinContent(e);
1418   Double_t v6 = fDiffFlow[rp][1][2]->GetBinContent(e);
1419   Double_t v8 = fDiffFlow[rp][1][3]->GetBinContent(e);
1420   // Errors:
1421   Double_t v2Error = fDiffFlow[rp][1][0]->GetBinError(e);
1422   Double_t v4Error = fDiffFlow[rp][1][1]->GetBinError(e);
1423   Double_t v6Error = fDiffFlow[rp][1][2]->GetBinError(e);
1424   Double_t v8Error = fDiffFlow[rp][1][3]->GetBinError(e);
1425   // Fill common hist results:
1426   if(rpPoi == "RP")
1427   {
1428    fCommonHistsResults2nd->FillDifferentialFlowEtaRP(e,v2,v2Error);
1429    fCommonHistsResults4th->FillDifferentialFlowEtaRP(e,v4,v4Error);
1430    fCommonHistsResults6th->FillDifferentialFlowEtaRP(e,v6,v6Error);
1431    fCommonHistsResults8th->FillDifferentialFlowEtaRP(e,v8,v8Error);
1432   } else if(rpPoi == "POI")
1433     {
1434      fCommonHistsResults2nd->FillDifferentialFlowEtaPOI(e,v2,v2Error);
1435      fCommonHistsResults4th->FillDifferentialFlowEtaPOI(e,v4,v4Error);
1436      fCommonHistsResults6th->FillDifferentialFlowEtaPOI(e,v6,v6Error);
1437      fCommonHistsResults8th->FillDifferentialFlowEtaPOI(e,v8,v8Error);
1438     }
1439  } // end of for(Int_t e=1;e<=fnBinsEta;e++)
1440
1441 } // end of void AliFlowAnalysisWithCumulants::FillCommonHistResultsForDifferentialFlow(TString rpPoi)
1442
1443 //================================================================================================================
1444
1445 void AliFlowAnalysisWithCumulants::CalculateDifferentialFlow(TString rpPoi, TString ptEta)
1446 {
1447  // Calculate differential flow for RPs/POIs vs pt/eta from cumulants.
1448
1449  Int_t rp = -1; // RP or POI
1450  Int_t pe = -1; // pt or eta
1451
1452  if(rpPoi == "RP")
1453  {
1454   rp = 0;
1455  } else if(rpPoi == "POI")
1456    {
1457     rp = 1;
1458    }
1459  if(ptEta == "pt")
1460  {
1461   pe = 0;
1462  } else if(ptEta == "eta")
1463    {
1464     pe = 1;
1465    }
1466
1467  // Reference cumulants:
1468  Double_t gfc2 = fReferenceFlowCumulants->GetBinContent(1); // reference 2nd order cumulant
1469  Double_t gfc4 = fReferenceFlowCumulants->GetBinContent(2); // reference 4th order cumulant
1470  Double_t gfc6 = fReferenceFlowCumulants->GetBinContent(3); // reference 6th order cumulant
1471  Double_t gfc8 = fReferenceFlowCumulants->GetBinContent(4); // reference 8th order cumulant
1472
1473  Int_t nBins = fDiffFlowCumulants[rp][pe][0]->GetXaxis()->GetNbins();
1474
1475  for(Int_t b=1;b<=nBins;b++)
1476  {
1477   // Differential cumulants:
1478   Double_t gfd2 = fDiffFlowCumulants[rp][pe][0]->GetBinContent(b); // differential 2nd order cumulant
1479   Double_t gfd4 = fDiffFlowCumulants[rp][pe][1]->GetBinContent(b); // differential 4th order cumulant
1480   Double_t gfd6 = fDiffFlowCumulants[rp][pe][2]->GetBinContent(b); // differential 6th order cumulant
1481   Double_t gfd8 = fDiffFlowCumulants[rp][pe][3]->GetBinContent(b); // differential 8th order cumulant
1482   // Differential flow:
1483   Double_t v2 = 0.; // v'{2,GFC}
1484   Double_t v4 = 0.; // v'{4,GFC}
1485   Double_t v6 = 0.; // v'{6,GFC}
1486   Double_t v8 = 0.; // v'{8,GFC}
1487   // 2nd order:
1488   if(gfc2>0.)
1489   {
1490    v2 = gfd2/pow(gfc2,0.5);
1491    fDiffFlow[rp][pe][0]->SetBinContent(b,v2);
1492   }
1493   // 4th order:
1494   if(gfc4<0.)
1495   {
1496    v4 = -gfd4/pow(-gfc4,.75);
1497    fDiffFlow[rp][pe][1]->SetBinContent(b,v4);
1498   }
1499   // 6th order:
1500   if(gfc6>0.)
1501   {
1502    v6 = gfd6/(4.*pow((1./4.)*gfc6,(5./6.)));
1503    fDiffFlow[rp][pe][2]->SetBinContent(b,v6);
1504   }
1505   // 8th order:
1506   if(gfc8<0.)
1507   {
1508    v8 = -gfd8/(33.*pow(-(1./33.)*gfc8,(7./8.)));
1509    fDiffFlow[rp][pe][3]->SetBinContent(b,v8);
1510   }
1511  } // end of for(Int_t b=1;b<=nBins;b++)
1512
1513 } // end of void AliFlowAnalysisWithCumulants::CalculateDifferentialFlow(TString rpPoi,TString ptEta)
1514
1515 //================================================================================================================
1516
1517 void AliFlowAnalysisWithCumulants::CalculateDifferentialFlowErrors(TString rpPoi,TString ptEta)
1518 {
1519  // Calculate errors of differential flow.
1520
1521  Int_t rp = -1; // RP or POI
1522  Int_t pe = -1; // pt or eta
1523
1524  if(rpPoi == "RP")
1525  {
1526   rp = 0;
1527  } else if(rpPoi == "POI")
1528    {
1529     rp = 1;
1530    }
1531  if(ptEta == "pt")
1532  {
1533   pe = 0;
1534  } else if(ptEta == "eta")
1535    {
1536     pe = 1;
1537    }
1538
1539  // Resolution chi:
1540  Double_t chi2 = fChi->GetBinContent(1);
1541  Double_t chi4 = fChi->GetBinContent(2);
1542  //Double_t chi6 = fChi->GetBinContent(3);
1543  //Double_t chi8 = fChi->GetBinContent(4);
1544
1545  Int_t nBins = fNoOfParticlesInBin[rp][pe]->GetXaxis()->GetNbins();
1546  for(Int_t b=1;b<=nBins;b++)
1547  {
1548   Int_t nParticles = (Int_t)fNoOfParticlesInBin[rp][pe]->GetBinEntries(b);
1549   // Error of 2nd order estimate:
1550   if(chi2>0. &&  nParticles>0.)
1551   {
1552    Double_t v2Error = pow((1./(2.*nParticles))*((1.+pow(chi2,2.))/pow(chi2,2.)),0.5);
1553    fDiffFlow[rp][pe][0]->SetBinError(b,v2Error);
1554   }
1555   // Error of 4th order estimate:
1556   if(chi4>0. &&  nParticles>0.)
1557   {
1558    Double_t v4Error = pow((1./(2.*nParticles))*((2.+6.*pow(chi4,2.)+pow(chi4,4.)+pow(chi4,6.))/pow(chi4,6.)),0.5);
1559    fDiffFlow[rp][pe][1]->SetBinError(b,v4Error);
1560   }
1561   // Error of 6th order estimate:
1562   //if(chi6>0. &&  nParticles>0.)
1563   //{
1564   // Double_t v6Error = ... // to be improved - yet to be calculated
1565    fDiffFlow[rp][pe][2]->SetBinError(b,0.);
1566   //}
1567   // Error of 8th order estimate:
1568   //if(chi8>0. &&  nParticles>0.)
1569   //{
1570   // Double_t v8Error = ... // to be improved - yet to be calculated
1571    fDiffFlow[rp][pe][3]->SetBinError(b,0.);
1572   //}
1573  } // end of for(Int_t b=1;b<=nBins;b++)
1574
1575 } // end of void AliFlowAnalysisWithCumulants::CalculateDifferentialFlowErrors(TString rpPoi,TString ptEta)
1576
1577 //================================================================================================================
1578
1579 void AliFlowAnalysisWithCumulants::CalculateCumulantsForDiffFlow(TString rpPoi,TString ptEta)
1580 {
1581  // Calculate cumulants for differential flow.
1582
1583  Int_t rp = -1; // RP or POI
1584  Int_t pe = -1; // pt or eta
1585
1586  if(rpPoi == "RP")
1587  {
1588   rp = 0;
1589  } else if(rpPoi == "POI")
1590    {
1591     rp = 1;
1592    }
1593  if(ptEta == "pt")
1594  {
1595   pe = 0;
1596  } else if(ptEta == "eta")
1597    {
1598     pe = 1;
1599    }
1600
1601  // [nBins][pMax][qMax]:
1602  Int_t nBins = fDiffFlowGenFun[0][rp][pe]->GetXaxis()->GetNbins();
1603  Int_t pMax  = fDiffFlowGenFun[0][rp][pe]->GetYaxis()->GetNbins();
1604  Int_t qMax  = fDiffFlowGenFun[0][rp][pe]->GetZaxis()->GetNbins();
1605  // <G[p][q]>
1606  TMatrixD dAvG(pMax,qMax);
1607  dAvG.Zero();
1608  for(Int_t p=0;p<pMax;p++)
1609  {
1610   for(Int_t q=0;q<qMax;q++)
1611   {
1612    dAvG(p,q) = fReferenceFlowGenFun->GetBinContent(fReferenceFlowGenFun->GetBin(p+1,q+1));
1613   }
1614  }
1615  // Loop over pt/eta bins and calculate differential cumulants:
1616  for(Int_t b=0;b<nBins;b++)
1617  {
1618   Double_t gfc[5] = {0.}; // to be improved (hardwired 5)
1619   Double_t D[5] = {0.}; // D_{p} in Eq. (11) in Practical guide // to be improved (hardwired 5)
1620   // ptBinRPNoOfParticles[b]=fPtBinRPNoOfParticles->GetBinEntries(b+1);
1621   for(Int_t p=0;p<pMax;p++)
1622   {
1623    Double_t tempSum = 0.;
1624    for(Int_t q=0;q<qMax;q++)
1625    {
1626     if(TMath::Abs(dAvG(p,q))>1.e-44)
1627     {
1628      Double_t X = fDiffFlowGenFun[0][rp][pe]->GetBinContent(fDiffFlowGenFun[0][rp][pe]->GetBin(b+1,p+1,q+1))/dAvG(p,q); // see Ollitrault's Practical guide (Eq. 11)
1629      Double_t Y = fDiffFlowGenFun[1][rp][pe]->GetBinContent(fDiffFlowGenFun[0][rp][pe]->GetBin(b+1,p+1,q+1))/dAvG(p,q); // see Ollitrault's Practical guide (Eq. 11)
1630      tempSum += cos(fMultiple*2.*q*TMath::Pi()/qMax)*X
1631               + sin(fMultiple*2.*q*TMath::Pi()/qMax)*Y;
1632     }
1633    }
1634    D[p] = (pow(fR0*pow(p+1.0,0.5),fMultiple)/qMax)*tempSum;
1635   }
1636   gfc[0] = (1./(fR0*fR0))*(5.*D[0]-5.*D[1]+(10./3.)*D[2]-(5./4.)*D[3]+(1./5.)*D[4]);
1637   gfc[1] = (1./pow(fR0,4.))*((-77./6.)*D[0]+(107./6.)*D[1]-(13./1.)*D[2]+(61./12.)*D[3]-(5./6.)*D[4]);
1638   gfc[2] = (1./pow(fR0,6.))*((71./2.)*D[0]-59.*D[1]+49.*D[2]-(41./2.)*D[3]+(7./2.)*D[4]);
1639   gfc[3] = (1./pow(fR0,8.))*(-84.*D[0]+156.*D[1]-144.*D[2]+66.*D[3]-12.*D[4]);
1640   // gfc[4] = (1./pow(fR0,10.))*(120.*D[0]-240.*D[1]+240.*D[2]-120.*D[3]+24.*D[4]); // 10th order cumulant (to be improved - where to store it?)
1641   // Store cumulants:
1642   for(Int_t co=0;co<4;co++)
1643   {
1644    fDiffFlowCumulants[rp][pe][co]->SetBinContent(b+1,gfc[co]);
1645   }
1646  }
1647
1648 } // end of void AliFlowAnalysisWithCumulants::CalculateCumulantsForDiffFlow(TString rpPoi, TString ptEta)
1649
1650 //================================================================================================================
1651
1652 void AliFlowAnalysisWithCumulants::PrintFinalResults(TString type)
1653 {
1654  // Printing on the screen the final results for reference flow and for integrated flow of RPs and POIs.
1655
1656  Int_t n = fHarmonic;
1657
1658  Double_t dVn[4] = {0.}; // array to hold Vn{2}, Vn{4}, Vn{6} and Vn{8}
1659  Double_t dVnErr[4] = {0.}; // array to hold errors of Vn{2}, Vn{4}, Vn{6} and Vn{8}
1660
1661  if(type == "RF")
1662  {
1663   dVn[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinContent(1);
1664   dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinError(1);
1665   dVn[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinContent(1);
1666   dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinError(1);
1667   dVn[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinContent(1);
1668   dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinError(1);
1669   dVn[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinContent(1);
1670   dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinError(1);
1671  } else if(type == "RP")
1672    {
1673     dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinContent(1);
1674     dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinError(1);
1675     dVn[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinContent(1);
1676     dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinError(1);
1677     dVn[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinContent(1);
1678     dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinError(1);
1679     dVn[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinContent(1);
1680     dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinError(1);
1681    } else if(type == "POI")
1682      {
1683       dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinContent(1);
1684       dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinError(1);
1685       dVn[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinContent(1);
1686       dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinError(1);
1687       dVn[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinContent(1);
1688       dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinError(1);
1689       dVn[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinContent(1);
1690       dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinError(1);
1691      } else
1692        {
1693         cout<<endl;
1694         cout<<" WARNING: Impossible type (can be RF, RP or POI) !!!!"<<endl;
1695         cout<<"          Results will not be printed on the screen."<<endl;
1696         cout<<endl;
1697         exit(0);
1698        }
1699
1700  TString title = " flow estimates from GF-cumulants";
1701  TString subtitle = "      (";
1702
1703  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
1704  {
1705   subtitle.Append(type);
1706   subtitle.Append(", without weights)");
1707  } else
1708    {
1709     subtitle.Append(type);
1710     subtitle.Append(", with weights)");
1711    }
1712
1713  cout<<endl;
1714  cout<<"*************************************"<<endl;
1715  cout<<"*************************************"<<endl;
1716  cout<<title.Data()<<endl;
1717  cout<<subtitle.Data()<<endl;
1718  cout<<endl;
1719
1720  for(Int_t i=0;i<4;i++)
1721  {
1722   cout<<"  v_"<<n<<"{"<<2*(i+1)<<"} = "<<dVn[i]<<" +/- "<<dVnErr[i]<<endl;
1723  }
1724
1725  cout<<endl;
1726  if(type == "RF")
1727  {
1728   cout<<"     nEvts = "<<(Int_t)fCommonHists->GetHistMultRP()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()<<endl;
1729  }
1730  else if (type == "RP")
1731  {
1732   cout<<"     nEvts = "<<(Int_t)fCommonHists->GetHistMultRP()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()<<endl;
1733  }
1734  else if (type == "POI")
1735  {
1736   cout<<"     nEvts = "<<(Int_t)fCommonHists->GetHistMultPOI()->GetEntries()<<", <M> = "<<(Double_t)fCommonHists->GetHistMultPOI()->GetMean()<<endl;
1737  }
1738  cout<<"*************************************"<<endl;
1739  cout<<"*************************************"<<endl;
1740  cout<<endl;
1741
1742 } // end of AliFlowAnalysisWithCumulants::PrintFinalResults(TString type);
1743
1744 //================================================================================================================
1745
1746 void AliFlowAnalysisWithCumulants::FillCommonHistResultsForReferenceFlow()
1747 {
1748  // Fill in AliFlowCommonHistResults dedicated histograms for reference flow.
1749
1750  // Results:
1751  Double_t v2 = fReferenceFlow->GetBinContent(1);
1752  Double_t v4 = fReferenceFlow->GetBinContent(2);
1753  Double_t v6 = fReferenceFlow->GetBinContent(3);
1754  Double_t v8 = fReferenceFlow->GetBinContent(4);
1755  // Errors:
1756  Double_t v2Error = fReferenceFlow->GetBinError(1);
1757  Double_t v4Error = fReferenceFlow->GetBinError(2);
1758  Double_t v6Error = fReferenceFlow->GetBinError(3);
1759  Double_t v8Error = fReferenceFlow->GetBinError(4);
1760  // Fill results end errors in common hist results:
1761  fCommonHistsResults2nd->FillIntegratedFlow(v2,v2Error);
1762  fCommonHistsResults4th->FillIntegratedFlow(v4,v4Error);
1763  fCommonHistsResults6th->FillIntegratedFlow(v6,v6Error);
1764  fCommonHistsResults8th->FillIntegratedFlow(v8,v8Error);
1765  // Chi:
1766  Double_t chi2 = fChi->GetBinContent(1);
1767  Double_t chi4 = fChi->GetBinContent(2);
1768  Double_t chi6 = fChi->GetBinContent(3);
1769  Double_t chi8 = fChi->GetBinContent(4);
1770  // Fill resolution chi in common hist results:
1771  fCommonHistsResults2nd->FillChi(chi2);
1772  fCommonHistsResults4th->FillChi(chi4);
1773  fCommonHistsResults6th->FillChi(chi6);
1774  fCommonHistsResults8th->FillChi(chi8);
1775
1776 } // end of AliFlowAnalysisWithCumulants::FillCommonHistResultsForReferenceFlow()
1777
1778 //================================================================================================================
1779
1780 void AliFlowAnalysisWithCumulants::CalculateReferenceFlowError()
1781 {
1782  // Calculate error of reference flow harmonics.
1783
1784  // Generating Function Cumulants:
1785  Double_t gfc2 = fReferenceFlowCumulants->GetBinContent(1); // GFC{2}
1786  Double_t gfc4 = fReferenceFlowCumulants->GetBinContent(2); // GFC{4}
1787  Double_t gfc6 = fReferenceFlowCumulants->GetBinContent(3); // GFC{6}
1788  Double_t gfc8 = fReferenceFlowCumulants->GetBinContent(4); // GFC{8}
1789  // Reference flow estimates:
1790  Double_t v2 = fReferenceFlow->GetBinContent(1); // v{2,GFC}
1791  Double_t v4 = fReferenceFlow->GetBinContent(2); // v{4,GFC}
1792  Double_t v6 = fReferenceFlow->GetBinContent(3); // v{6,GFC}
1793  Double_t v8 = fReferenceFlow->GetBinContent(4); // v{8,GFC}
1794  // Statistical errors of reference flow estimates:
1795  Double_t v2Error = 0.; // statistical error of v{2,GFC}
1796  Double_t v4Error = 0.; // statistical error of v{4,GFC}
1797  Double_t v6Error = 0.; // statistical error of v{6,GFC}
1798  Double_t v8Error = 0.; // statistical error of v{8,GFC}
1799  // Chi:
1800  Double_t chi2 = 0.;
1801  Double_t chi4 = 0.;
1802  Double_t chi6 = 0.;
1803  Double_t chi8 = 0.;
1804  // <Q-vector stuff>:
1805  Double_t dAvQx  = fQvectorComponents->GetBinContent(1); // <Q_x>
1806  Double_t dAvQy  = fQvectorComponents->GetBinContent(2); // <Q_y>
1807  Double_t dAvQ2x = fQvectorComponents->GetBinContent(3); // <(Q_x)^2>
1808  Double_t dAvQ2y = fQvectorComponents->GetBinContent(4); // <(Q_y)^2>
1809  // <w^2>:
1810  Double_t dAvw2 = 1.;
1811  if(fnEvts>0)
1812  {
1813   dAvw2 = fAverageOfSquaredWeight->GetBinContent(1);
1814   if(TMath::Abs(dAvw2)<1.e-44)
1815   {
1816    cout<<endl;
1817    cout<<" WARNING (GFC): Average of squared weight is 0 in GFC. Most probably one of the histograms"<<endl;
1818    cout<<"                in the file \"weights.root\" was empty. Nothing will be calculated !!!!"<<endl;
1819    cout<<endl;
1820   }
1821  }
1822  // Calculating statistical error of v{2,GFC}:
1823  if(fnEvts>0. && fAvM>0. && dAvw2>0. && gfc2>=0.)
1824  {
1825   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(gfc2,(1./2.))*(fAvM/dAvw2),2.)>0.))
1826   {
1827    chi2 = (fAvM*v2)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v2*fAvM/dAvw2,2.),0.5);
1828   }
1829   if(TMath::Abs(chi2)>1.e-44)
1830   {
1831    v2Error = pow(((1./(2.*fAvM*fnEvts))*((1.+2.*pow(chi2,2))/(2.*pow(chi2,2)))),0.5);
1832   }
1833  }
1834  // Calculating statistical error of v{4,GFC}:
1835  if(fnEvts>0 && fAvM>0 && dAvw2>0 && gfc4<=0.)
1836  {
1837   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-gfc4,(1./4.))*(fAvM/dAvw2),2.)>0.))
1838   {
1839    chi4 = (fAvM*v4)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v4*fAvM/dAvw2,2.),0.5);
1840   }
1841   if(TMath::Abs(chi4)>1.e-44)
1842   {
1843    v4Error = (1./(pow(2.*fAvM*fnEvts,0.5)))*pow((1.+4.*pow(chi4,2)+1.*pow(chi4,4.)+2.*pow(chi4,6.))/(2.*pow(chi4,6.)),0.5);
1844   }
1845  }
1846  // Calculating statistical error of v{6,GFC}:
1847  if(fnEvts>0 && fAvM>0 && dAvw2>0 && gfc6>=0.)
1848  {
1849   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*gfc6,(1./6.))*(fAvM/dAvw2),2.)>0.))
1850   {
1851    chi6 = (fAvM*v6)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v6*fAvM/dAvw2,2.),0.5);
1852   }
1853   if(TMath::Abs(chi6)>1.e-44)
1854   {
1855    v6Error = (1./(pow(2.*fAvM*fnEvts,0.5)))*pow((3.+18.*pow(chi6,2)+9.*pow(chi6,4.)+28.*pow(chi6,6.)
1856               +12.*pow(chi6,8.)+24.*pow(chi6,10.))/(24.*pow(chi6,10.)),0.5);
1857   }
1858  }
1859  // Calculating statistical error of v{8,GFC}:
1860  if(fnEvts>0 && fAvM>0 && dAvw2>0 && gfc8<=0.)
1861  {
1862   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-(1./33.)*gfc8,(1./8.))*(fAvM/dAvw2),2.)>0.))
1863   {
1864    chi8=(fAvM*v8)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v8*fAvM/dAvw2,2.),0.5);
1865   }
1866   if(TMath::Abs(chi8)>1.e-44)
1867   {
1868    v8Error = (1./(pow(2.*fAvM*fnEvts,0.5)))*pow((12.+96.*pow(chi8,2.)+72.*pow(chi8,4.)+304.*pow(chi8,6.)
1869               +257.*pow(chi8,8.)+804.*pow(chi8,10.)+363.*pow(chi8,12.)+726.*pow(chi8,14.))/(726.*pow(chi8,14.)),0.5);
1870   }
1871  }
1872
1873  // Store errors for reference flow:
1874  fReferenceFlow->SetBinError(1,v2Error);
1875  fReferenceFlow->SetBinError(2,v4Error);
1876  fReferenceFlow->SetBinError(3,v6Error);
1877  fReferenceFlow->SetBinError(4,v8Error);
1878  // Store resolution chi:
1879  fChi->SetBinContent(1,chi2);
1880  fChi->SetBinContent(2,chi4);
1881  fChi->SetBinContent(3,chi6);
1882  fChi->SetBinContent(4,chi8);
1883
1884 } // end of void AliFlowAnalysisWithCumulants::CalculateReferenceFlowError()
1885
1886 //================================================================================================================
1887
1888 void AliFlowAnalysisWithCumulants::CalculateReferenceFlow()
1889 {
1890  // Calculate from isotropic cumulants reference flow.
1891
1892  // Generating Function Cumulants:
1893  Double_t gfc2 = fReferenceFlowCumulants->GetBinContent(1); // GFC{2}
1894  Double_t gfc4 = fReferenceFlowCumulants->GetBinContent(2); // GFC{4}
1895  Double_t gfc6 = fReferenceFlowCumulants->GetBinContent(3); // GFC{6}
1896  Double_t gfc8 = fReferenceFlowCumulants->GetBinContent(4); // GFC{8}
1897  // Double_t gfc10 = fReferenceFlowCumulants->GetBinContent(5); // GFC{10} overflow
1898  // Reference flow estimates:
1899  Double_t v2 = 0.; // v{2,GFC}
1900  Double_t v4 = 0.; // v{4,GFC}
1901  Double_t v6 = 0.; // v{6,GFC}
1902  Double_t v8 = 0.; // v{8,GFC}
1903  // Double_t v10 = 0.; // v{10,GFC} overflow
1904  // Calculate reference flow estimates from Q-cumulants:
1905  if(gfc2>=0.) v2 = pow(gfc2,1./2.);
1906  if(gfc4<=0.) v4 = pow(-1.*gfc4,1./4.);
1907  if(gfc6>=0.) v6 = pow((1./4.)*gfc6,1./6.);
1908  if(gfc8<=0.) v8 = pow((-1./33.)*gfc8,1./8.);
1909  // if(gfc10>=0.) v10 = pow((1./456.)*gfc10,1./10.);
1910  // Store results for reference flow:
1911  fReferenceFlow->SetBinContent(1,v2);
1912  fReferenceFlow->SetBinContent(2,v4);
1913  fReferenceFlow->SetBinContent(3,v6);
1914  fReferenceFlow->SetBinContent(4,v8);
1915  // fReferenceFlow->SetBinContent(5,v10); // overflow
1916
1917 } // end of void AliFlowAnalysisWithCumulants::CalculateReferenceFlow()
1918
1919 //================================================================================================================
1920
1921 void AliFlowAnalysisWithCumulants::CalculateCumulantsForReferenceFlow()
1922 {
1923  // Calculate cumulants for reference flow.
1924
1925  Int_t pMax = fReferenceFlowGenFun->GetXaxis()->GetNbins();
1926  Int_t qMax = fReferenceFlowGenFun->GetYaxis()->GetNbins();
1927
1928  // <G[p][q]>
1929  TMatrixD dAvG(pMax,qMax);
1930  dAvG.Zero();
1931  Bool_t someAvGEntryIsNegative = kFALSE;
1932  for(Int_t p=0;p<pMax;p++)
1933  {
1934   for(Int_t q=0;q<qMax;q++)
1935   {
1936    dAvG(p,q) = fReferenceFlowGenFun->GetBinContent(fReferenceFlowGenFun->GetBin(p+1,q+1));
1937    if(dAvG(p,q)<0.)
1938    {
1939     someAvGEntryIsNegative = kTRUE;
1940     cout<<endl;
1941     cout<<" WARNING: "<<Form("<G[%d][%d]> is negative !!!! GFC results are meaningless.",p,q)<<endl;
1942     cout<<endl;
1943    }
1944   }
1945  }
1946
1947  // C[p][q] (generating function for the cumulants)
1948  TMatrixD dC(pMax,qMax);
1949  dC.Zero();
1950  if(fAvM>0. && !someAvGEntryIsNegative)
1951  {
1952   for(Int_t p=0;p<pMax;p++)
1953   {
1954    for(Int_t q=0;q<qMax;q++)
1955    {
1956     dC(p,q) = fAvM*(pow(dAvG(p,q),(1./fAvM))-1.);
1957    }
1958   }
1959  }
1960
1961  // Averaging the generating function for cumulants over azimuth
1962  // in order to eliminate detector effects.
1963  // <C[p][q]> (Remark: here <> stands for average over azimuth):
1964  TVectorD dAvC(pMax);
1965  dAvC.Zero();
1966  for(Int_t p=0;p<pMax;p++)
1967  {
1968   Double_t temp = 0.;
1969   for(Int_t q=0;q<qMax;q++)
1970   {
1971    temp += 1.*dC(p,q);
1972   }
1973   dAvC[p] = temp/qMax;
1974  }
1975
1976  // Finally, the isotropic cumulants for reference flow:
1977  TVectorD cumulant(pMax);
1978  cumulant.Zero();
1979  cumulant[0] = (-1./(60*fR0*fR0))*((-300.)*dAvC[0]+300.*dAvC[1]-200.*dAvC[2]+75.*dAvC[3]-12.*dAvC[4]);
1980  cumulant[1] = (-1./(6.*pow(fR0,4.)))*(154.*dAvC[0]-214.*dAvC[1]+156.*dAvC[2]-61.*dAvC[3]+10.*dAvC[4]);
1981  cumulant[2] = (3./(2.*pow(fR0,6.)))*(71.*dAvC[0]-118.*dAvC[1]+98.*dAvC[2]-41.*dAvC[3]+7.*dAvC[4]);
1982  cumulant[3] = (-24./pow(fR0,8.))*(14.*dAvC[0]-26.*dAvC[1]+24.*dAvC[2]-11.*dAvC[3]+2.*dAvC[4]);
1983  cumulant[4] = (120./pow(fR0,10.))*(5.*dAvC[0]-10.*dAvC[1]+10.*dAvC[2]-5.*dAvC[3]+1.*dAvC[4]);
1984
1985  // Store cumulants:
1986  // Remark: the highest order cumulant is on purpose in the overflow.
1987  for(Int_t co=0;co<pMax;co++) // cumulant order
1988  {
1989   fReferenceFlowCumulants->SetBinContent(co+1,cumulant[co]);
1990  }
1991
1992 } // end of void AliFlowAnalysisWithCumulants::CalculateCumulantsForReferenceFlow()
1993
1994 //================================================================================================================
1995
1996 void AliFlowAnalysisWithCumulants::GetAvMultAndNoOfEvts()
1997 {
1998  // From relevant common control histogram get average multiplicity of RPs and number of events.
1999
2000  fAvM = (Double_t)fCommonHists->GetHistMultRP()->GetMean();
2001  fnEvts = (Int_t)fCommonHists->GetHistMultRP()->GetEntries();
2002
2003 } // end of void AliFlowAnalysisWithCumulants::GetAvMultAndNoOfEvts()
2004
2005 //================================================================================================================
2006
2007 void AliFlowAnalysisWithCumulants::InitializeArrays()
2008 {
2009  // Initialize all arrays.
2010
2011  for(Int_t ri=0;ri<2;ri++)
2012  {
2013   for(Int_t rp=0;rp<2;rp++)
2014   {
2015    for(Int_t pe=0;pe<2;pe++)
2016    {
2017     fDiffFlowGenFun[ri][rp][pe] = NULL;
2018    }
2019   }
2020  }
2021  for(Int_t rp=0;rp<2;rp++)
2022  {
2023   for(Int_t pe=0;pe<2;pe++)
2024   {
2025    fNoOfParticlesInBin[rp][pe] = NULL;
2026   }
2027  }
2028  for(Int_t rp=0;rp<2;rp++)
2029  {
2030   for(Int_t pe=0;pe<2;pe++)
2031   {
2032    for(Int_t co=0;co<4;co++)
2033    {
2034     fDiffFlowCumulants[rp][pe][co] = NULL;
2035     fDiffFlow[rp][pe][co] = NULL;
2036    }
2037   }
2038  }
2039  for(Int_t i=0;i<3;i++)
2040  {
2041   fPrintFinalResults[i] = kTRUE;
2042  }
2043  for(Int_t r=0;r<10;r++)
2044  {
2045   fTuningR0[r] = 0.;
2046   for(Int_t pq=0;pq<5;pq++)
2047   {
2048    fTuningGenFun[r][pq] = NULL;
2049    fTuningCumulants[r][pq] = NULL;
2050    fTuningFlow[r][pq] = NULL;
2051   }
2052  }
2053
2054 } // end of void AliFlowAnalysisWithCumulants::InitializeArrays()
2055
2056 //================================================================================================================
2057
2058 void AliFlowAnalysisWithCumulants::CrossCheckSettings()
2059 {
2060  // Cross-check the user settings before starting.
2061
2062  // a) Cross check if the choice for multiplicity weight make sense.
2063
2064  // a) Cross check if the choice for multiplicity weight make sense:
2065  if(strcmp(fMultiplicityWeight->Data(),"unit") &&
2066     strcmp(fMultiplicityWeight->Data(),"multiplicity"))
2067  {
2068   cout<<endl;
2069   cout<<"WARNING (GFC): Multiplicity weight can be either \"unit\" or \"multiplicity\"."<<endl;
2070   cout<<"               Certainly not \""<<fMultiplicityWeight->Data()<<"\"."<<endl;
2071   cout<<endl;
2072   exit(0);
2073  }
2074
2075
2076
2077 } // end of void AliFlowAnalysisWithCumulants::CrossCheckSettings()
2078
2079 //================================================================================================================
2080
2081 void AliFlowAnalysisWithCumulants::AccessConstants()
2082 {
2083  // Access needed common constants from AliFlowCommonConstants.
2084
2085  fnBinsPhi = AliFlowCommonConstants::GetMaster()->GetNbinsPhi();
2086  fPhiMin = AliFlowCommonConstants::GetMaster()->GetPhiMin();
2087  fPhiMax = AliFlowCommonConstants::GetMaster()->GetPhiMax();
2088  if(fnBinsPhi) fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi;
2089  fnBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
2090  fPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
2091  fPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
2092  if(fnBinsPt) fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt;
2093  fnBinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
2094  fEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
2095  fEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
2096  if(fnBinsEta) fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta;
2097
2098 } // end of void AliFlowAnalysisWithCumulants::AccessConstants()
2099
2100 //================================================================================================================
2101
2102 void AliFlowAnalysisWithCumulants::BookAndFillWeightsHistograms()
2103 {
2104  // book and fill histograms which hold phi, pt and eta weights
2105
2106  if(!fWeightsList)
2107  {
2108   cout<<"WARNING (GFC): fWeightsList is NULL in AFAWGFC::BAFWH() !!!!"<<endl;
2109   exit(0);
2110  }
2111
2112  TString fUseParticleWeightsName = "fUseParticleWeights";
2113  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
2114  fUseParticleWeights->SetLabelSize(0.06);
2115  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
2116  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
2117  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
2118  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
2119  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
2120  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
2121  fWeightsList->Add(fUseParticleWeights);
2122  if(fUsePhiWeights)
2123  {
2124   if(fWeightsList->FindObject("phi_weights"))
2125   {
2126    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
2127    if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.))
2128    {
2129     cout<<endl;
2130     cout<<"WARNING (GFC): Inconsistent binning in histograms for phi-weights throughout the code."<<endl;
2131     cout<<endl;
2132     exit(0);
2133    }
2134   } else
2135     {
2136      cout<<endl;
2137      cout<<"WARNING (GFC): fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWGFC::BAFWH() !!!!"<<endl;
2138      cout<<endl;
2139      exit(0);
2140     }
2141  } // end of if(fUsePhiWeights)
2142
2143  if(fUsePtWeights)
2144  {
2145   if(fWeightsList->FindObject("pt_weights"))
2146   {
2147    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
2148    if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.))
2149    {
2150     cout<<endl;
2151     cout<<"WARNING (GFC): Inconsistent binning in histograms for pt-weights throughout the code."<<endl;
2152     cout<<endl;
2153     exit(0);
2154    }
2155   } else
2156     {
2157      cout<<endl;
2158      cout<<"WARNING (GFC): fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWGFC::BAFWH() !!!!"<<endl;
2159      cout<<endl;
2160      exit(0);
2161     }
2162  } // end of if(fUsePtWeights)
2163
2164  if(fUseEtaWeights)
2165  {
2166   if(fWeightsList->FindObject("eta_weights"))
2167   {
2168    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
2169    if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.))
2170    {
2171     cout<<endl;
2172     cout<<"WARNING (GFC): Inconsistent binning in histograms for eta-weights throughout the code."<<endl;
2173     cout<<endl;
2174     exit(0);
2175    }
2176   } else
2177     {
2178      cout<<endl;
2179      cout<<"WARNING (GFC): fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWGFC::BAFWH() !!!!"<<endl;
2180      cout<<endl;
2181      exit(0);
2182     }
2183  } // end of if(fUseEtaWeights)
2184
2185 } // end of AliFlowAnalysisWithCumulants::BookAndFillWeightsHistograms()
2186
2187 //================================================================================================================
2188
2189 void AliFlowAnalysisWithCumulants::BookEverythingForReferenceFlow()
2190 {
2191  // Book all objects relevant for calculation of reference flow.
2192
2193  // a) Define static constants for array's boundaries;
2194  // b) Book profile to hold all flags for reference flow;
2195  // c) Book all event-by-event quantities;
2196  // d) Book all profiles;
2197  // e) Book all histograms.
2198
2199  // a) Define static constants for array's boundaries:
2200  static const Int_t pMax = 5;
2201  static const Int_t qMax = 11;
2202
2203  // b) Book profile to hold all flags for reference flow:
2204  TString referenceFlowFlagsName = "fReferenceFlowFlags";
2205  fReferenceFlowFlags = new TProfile(referenceFlowFlagsName.Data(),"Flags for Reference Flow",2,0,2);
2206  fReferenceFlowFlags->SetTickLength(-0.01,"Y");
2207  fReferenceFlowFlags->SetMarkerStyle(25);
2208  fReferenceFlowFlags->SetLabelSize(0.05);
2209  fReferenceFlowFlags->SetLabelOffset(0.02,"Y");
2210  fReferenceFlowFlags->GetXaxis()->SetBinLabel(1,"Particle weights");
2211  fReferenceFlowFlags->GetXaxis()->SetBinLabel(2,"Event weights");
2212  fReferenceFlowList->Add(fReferenceFlowFlags);
2213
2214  // c) Book all event-by-event quantities:
2215  fGEBE = new TMatrixD(pMax,qMax);
2216
2217  // d) Book all profiles:
2218  // Average of the generating function for reference flow <G[p][q]>:
2219  fReferenceFlowGenFun = new TProfile2D("fReferenceFlowGenFun","#LTG[p][q]#GT",pMax,0.,(Double_t)pMax,qMax,0.,(Double_t)qMax);
2220  fReferenceFlowGenFun->SetXTitle("p");
2221  fReferenceFlowGenFun->SetYTitle("q");
2222  fReferenceFlowProfiles->Add(fReferenceFlowGenFun);
2223  // Averages of Q-vector components:
2224  fQvectorComponents = new TProfile("fQvectorComponents","Averages of Q-vector components",4,0.,4.);
2225  fQvectorComponents->SetLabelSize(0.06);
2226  fQvectorComponents->SetMarkerStyle(25);
2227  fQvectorComponents->GetXaxis()->SetBinLabel(1,"#LTQ_{x}#GT"); // Q_{x}
2228  fQvectorComponents->GetXaxis()->SetBinLabel(2,"#LTQ_{y}#GT"); // Q_{y}
2229  fQvectorComponents->GetXaxis()->SetBinLabel(3,"#LTQ_{x}^{2}#GT"); // Q_{x}^{2}
2230  fQvectorComponents->GetXaxis()->SetBinLabel(4,"#LTQ_{y}^{2}#GT"); // Q_{y}^{2}
2231  fReferenceFlowProfiles->Add(fQvectorComponents);
2232  // <<w^2>>, where w = wPhi*wPt*wEta:
2233  fAverageOfSquaredWeight = new TProfile("fAverageOfSquaredWeight","#LT#LTw^{2}#GT#GT",1,0,1);
2234  fAverageOfSquaredWeight->SetLabelSize(0.06);
2235  fAverageOfSquaredWeight->SetMarkerStyle(25);
2236  fAverageOfSquaredWeight->SetLabelOffset(0.01);
2237  fAverageOfSquaredWeight->GetXaxis()->SetBinLabel(1,"#LT#LTw^{2}#GT#GT");
2238  fReferenceFlowProfiles->Add(fAverageOfSquaredWeight);
2239
2240  // e) Book all histograms:
2241  // Final results for isotropic cumulants for reference flow:
2242  TString referenceFlowCumulantsName = "fReferenceFlowCumulants";
2243  fReferenceFlowCumulants = new TH1D(referenceFlowCumulantsName.Data(),"Isotropic Generating Function Cumulants for reference flow",4,0,4); // to be improved (hw 4)
2244  fReferenceFlowCumulants->SetLabelSize(0.05);
2245  fReferenceFlowCumulants->SetMarkerStyle(25);
2246  fReferenceFlowCumulants->GetXaxis()->SetBinLabel(1,"GFC{2}");
2247  fReferenceFlowCumulants->GetXaxis()->SetBinLabel(2,"GFC{4}");
2248  fReferenceFlowCumulants->GetXaxis()->SetBinLabel(3,"GFC{6}");
2249  fReferenceFlowCumulants->GetXaxis()->SetBinLabel(4,"GFC{8}");
2250  fReferenceFlowResults->Add(fReferenceFlowCumulants);
2251  // Final results for reference flow:
2252  fReferenceFlow = new TH1D("fReferenceFlow","Reference flow",4,0,4); // to be improved (hardwired 4)
2253  fReferenceFlow->SetLabelSize(0.06);
2254  fReferenceFlow->SetMarkerStyle(25);
2255  fReferenceFlow->GetXaxis()->SetBinLabel(1,"v_{n}{2,GFC}");
2256  fReferenceFlow->GetXaxis()->SetBinLabel(2,"v_{n}{4,GFC}");
2257  fReferenceFlow->GetXaxis()->SetBinLabel(3,"v_{n}{6,GFC}");
2258  fReferenceFlow->GetXaxis()->SetBinLabel(4,"v_{n}{8,GFC}");
2259  fReferenceFlowResults->Add(fReferenceFlow);
2260  // Final results for resolution:
2261  fChi = new TH1D("fChi","Resolution",4,0,4); // to be improved (hardwired 4)
2262  fChi->SetLabelSize(0.06);
2263  fChi->SetMarkerStyle(25);
2264  fChi->GetXaxis()->SetBinLabel(1,"#chi_{2}");
2265  fChi->GetXaxis()->SetBinLabel(2,"#chi_{4}");
2266  fChi->GetXaxis()->SetBinLabel(3,"#chi_{6}");
2267  fChi->GetXaxis()->SetBinLabel(4,"#chi_{8}");
2268  fReferenceFlowResults->Add(fChi);
2269
2270 } // end of void AliFlowAnalysisWithCumulants::BookEverythingForReferenceFlow()
2271
2272 //================================================================================================================
2273
2274 void AliFlowAnalysisWithCumulants::BookEverythingForTuning()
2275 {
2276  // Book all objects relevant for tuning.
2277
2278  // a) Define pMax's and qMax's:
2279  // b) Book profile to hold all tuning parameters and flags;
2280  // c) Book all profiles;
2281  // d) Book all histograms.
2282
2283  // a) Define pMax's and qMax's:
2284  Int_t pMax[5] = {2,3,4,5,8};
2285  Int_t qMax[5] = {5,7,9,11,17};
2286
2287  // b) Book profile to hold all tuning parameters and flags:
2288  TString tuningFlagsName = "fTuningFlags";
2289  fTuningFlags = new TProfile(tuningFlagsName.Data(),"Tuning parameters",10,0,10);
2290  // fTuningFlags->SetTickLength(-0.01,"Y");
2291  fTuningFlags->SetMarkerStyle(25);
2292  fTuningFlags->SetLabelSize(0.05);
2293  fTuningFlags->SetLabelOffset(0.02,"X");
2294  for(Int_t r=1;r<=10;r++)
2295  {
2296   fTuningFlags->GetXaxis()->SetBinLabel(r,Form("r_{0,%d}",r-1));
2297   fTuningFlags->Fill(r-0.5,fTuningR0[r-1],1.);
2298  }
2299  fTuningList->Add(fTuningFlags);
2300
2301  // c) Book all profiles:
2302  // Average of the generating function for reference flow <G[p][q]> for different tuning parameters:
2303  for(Int_t r=0;r<10;r++)
2304  {
2305   for(Int_t pq=0;pq<5;pq++)
2306   {
2307    fTuningGenFun[r][pq] = new TProfile2D(Form("fTuningGenFun (r_{0,%i}, pq set %i)",r,pq),
2308                                          Form("#LTG[p][q]#GT for r_{0} = %f, p_{max} = %i, q_{max} = %i",fTuningR0[r],pMax[pq],qMax[pq]),
2309                                          pMax[pq],0.,(Double_t)pMax[pq],qMax[pq],0.,(Double_t)qMax[pq]);
2310    fTuningGenFun[r][pq]->SetXTitle("p");
2311    fTuningGenFun[r][pq]->SetYTitle("q");
2312    fTuningProfiles->Add(fTuningGenFun[r][pq]);
2313   }
2314  }
2315
2316  // d) Book all histograms:
2317  // Final results for isotropic cumulants for reference flow for different tuning parameters:
2318  for(Int_t r=0;r<10;r++)
2319  {
2320   for(Int_t pq=0;pq<5;pq++)
2321   {
2322    fTuningCumulants[r][pq] = new TH1D(Form("fTuningCumulants (r_{0,%i}, pq set %i)",r,pq),
2323                                       Form("GFC for r_{0} = %f, p_{max} = %i, q_{max} = %i",fTuningR0[r],pMax[pq],qMax[pq]),
2324                                       pMax[pq],0,pMax[pq]);
2325    // fTuningCumulants[r][pq]->SetLabelSize(0.05);
2326    fTuningCumulants[r][pq]->SetMarkerStyle(25);
2327    for(Int_t b=1;b<=pMax[pq];b++)
2328    {
2329     fTuningCumulants[r][pq]->GetXaxis()->SetBinLabel(b,Form("GFC{%i}",2*b));
2330    }
2331    fTuningResults->Add(fTuningCumulants[r][pq]);
2332   }
2333  }
2334  // Final results for reference flow for different tuning parameters:
2335  for(Int_t r=0;r<10;r++)
2336  {
2337   for(Int_t pq=0;pq<5;pq++)
2338   {
2339    fTuningFlow[r][pq] = new TH1D(Form("fTuningFlow (r_{0,%i}, pq set %i)",r,pq),
2340                           Form("Reference flow for r_{0} = %f, p_{max} = %i, q_{max} = %i",fTuningR0[r],pMax[pq],qMax[pq]),
2341                           pMax[pq],0,pMax[pq]);
2342    // fTuningFlow[r][pq]->SetLabelSize(0.06);
2343    fTuningFlow[r][pq]->SetMarkerStyle(25);
2344    for(Int_t b=1;b<=pMax[pq];b++)
2345    {
2346     fTuningFlow[r][pq]->GetXaxis()->SetBinLabel(b,Form("v{%i,GFC}",2*b));
2347    }
2348    fTuningResults->Add(fTuningFlow[r][pq]);
2349   }
2350  }
2351
2352 } // end of void AliFlowAnalysisWithCumulants::BookEverythingForTuning()
2353
2354 //================================================================================================================
2355
2356 void AliFlowAnalysisWithCumulants::BookEverythingForDiffFlow()
2357 {
2358  // Book all objects relevant for calculation of differential flow.
2359
2360  // a) Define static constants for array's boundaries;
2361  // b) Define local variables and local flags for booking;
2362  // c) Book profile to hold all flags for differential flow;
2363  // d) Book all event-by-event quantities;
2364  // e) Book all profiles;
2365  // f) Book all histograms.
2366
2367  // a) Define static constants for array's boundaries:
2368  static const Int_t pMax = 5;
2369  static const Int_t qMax = 11;
2370
2371  // b) Define local variables and local flags for booking:
2372  Int_t nBinsPtEta[2] = {fnBinsPt,fnBinsEta};
2373  Double_t minPtEta[2] = {fPtMin,fEtaMin};
2374  Double_t maxPtEta[2] = {fPtMax,fEtaMax};
2375  TString reIm[2] = {"Re","Im"};
2376  TString rpPoi[2] = {"RP","POI"};
2377  TString ptEta[2] = {"p_{t}","#eta"};
2378  TString order[4] = {"2nd order","4th order","6th order","8th order"};
2379
2380  // c) Book profile to hold all flags for differential flow:
2381  TString diffFlowFlagsName = "fDiffFlowFlags";
2382  fDiffFlowFlags = new TProfile(diffFlowFlagsName.Data(),"Flags for Differential Flow",1,0,1);
2383  fDiffFlowFlags->SetTickLength(-0.01,"Y");
2384  fDiffFlowFlags->SetMarkerStyle(25);
2385  fDiffFlowFlags->SetLabelSize(0.05);
2386  fDiffFlowFlags->SetLabelOffset(0.02,"Y");
2387  fDiffFlowFlags->GetXaxis()->SetBinLabel(1,"...");
2388  fDiffFlowList->Add(fDiffFlowFlags);
2389
2390  // d) Book all event-by-event quantities:
2391  // ... (to be improved - perhaps not needed)
2392
2393  // e) Book all profiles:
2394  // Generating functions for differential flow:
2395  for(Int_t ri=0;ri<2;ri++)
2396  {
2397   for(Int_t rp=0;rp<2;rp++)
2398   {
2399    for(Int_t pe=0;pe<2;pe++)
2400    {
2401     fDiffFlowGenFun[ri][rp][pe] = new TProfile3D(Form("fDiffFlowGenFun (%s, %s, %s)",reIm[ri].Data(),rpPoi[rp].Data(),ptEta[pe].Data()),
2402                                                  Form("#LT%s[D[%s-bin][p][q]]#GT for %ss",reIm[ri].Data(),ptEta[pe].Data(),rpPoi[rp].Data()),
2403                                                  nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe],pMax,0.,(Double_t)pMax,qMax,0.,(Double_t)qMax);
2404     fDiffFlowGenFun[ri][rp][pe]->SetXTitle(ptEta[pe].Data());
2405     fDiffFlowGenFun[ri][rp][pe]->SetYTitle("p");
2406     fDiffFlowGenFun[ri][rp][pe]->SetZTitle("q");
2407     fDiffFlowGenFun[ri][rp][pe]->SetTitleOffset(1.44,"X");
2408     fDiffFlowGenFun[ri][rp][pe]->SetTitleOffset(1.44,"Y");
2409     fDiffFlowProfiles->Add(fDiffFlowGenFun[ri][rp][pe]);
2410     // to be improved - alternative // nBinsPtEta[pe],(Double_t)(fPtMin/fPtBinWidth),(Double_t)(fPtMax/fPtBinWidth),pMax,0.,(Double_t)pMax,qMax,0.,(Double_t)qMax);
2411    }
2412   }
2413  }
2414  // Number of particles in pt/eta bin for RPs/POIs:
2415  for(Int_t rp=0;rp<2;rp++)
2416  {
2417   for(Int_t pe=0;pe<2;pe++)
2418   {
2419    fNoOfParticlesInBin[rp][pe] = new TProfile(Form("fNoOfParticlesInBin (%s, %s)",rpPoi[rp].Data(),ptEta[pe].Data()),
2420                                               Form("Number of %ss per %s bin",rpPoi[rp].Data(),ptEta[pe].Data()),
2421                                               nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
2422    fNoOfParticlesInBin[rp][pe]->SetXTitle(ptEta[pe].Data());
2423    fDiffFlowProfiles->Add(fNoOfParticlesInBin[rp][pe]);
2424   }
2425  }
2426  // Differential cumulants per pt/eta bin for RPs/POIs:
2427  for(Int_t rp=0;rp<2;rp++)
2428  {
2429   for(Int_t pe=0;pe<2;pe++)
2430   {
2431    for(Int_t co=0;co<4;co++)
2432    {
2433     fDiffFlowCumulants[rp][pe][co] = new TH1D(Form("fDiffFlowCumulants (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data()),
2434                                               Form("Differential %s cumulant for %ss vs %s",order[co].Data(),rpPoi[rp].Data(),ptEta[pe].Data()),
2435                                               nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
2436     fDiffFlowCumulants[rp][pe][co]->SetXTitle(ptEta[pe].Data());
2437     fDiffFlowResults->Add(fDiffFlowCumulants[rp][pe][co]);
2438    }
2439   }
2440  }
2441  // Differential flow per pt/eta bin for RPs/POIs:
2442  for(Int_t rp=0;rp<2;rp++)
2443  {
2444   for(Int_t pe=0;pe<2;pe++)
2445   {
2446    for(Int_t co=0;co<4;co++)
2447    {
2448     fDiffFlow[rp][pe][co] = new TH1D(Form("fDiffFlow (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data()),
2449                                      Form("Differential flow from %s cumulant for %ss vs %s",order[co].Data(),rpPoi[rp].Data(),ptEta[pe].Data()),
2450                                      nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
2451     fDiffFlow[rp][pe][co]->SetXTitle(ptEta[pe].Data());
2452     fDiffFlowResults->Add(fDiffFlow[rp][pe][co]);
2453    }
2454   }
2455  }
2456
2457 }// end of void AliFlowAnalysisWithCumulants::BookEverythingForDiffFlow()
2458
2459 //================================================================================================================
2460
2461 void AliFlowAnalysisWithCumulants::StoreReferenceFlowFlags()
2462 {
2463  // Store all flags for reference flow in profile fReferenceFlowFlags.
2464
2465  if(!fReferenceFlowFlags)
2466  {
2467   cout<<endl;
2468   cout<<"WARNING: !fReferenceFlowFlags is NULL in AFAWC::SRFF() !!!!"<<endl;
2469   cout<<endl;
2470   exit(0);
2471  }
2472
2473  // Particle weights used or not:
2474  fReferenceFlowFlags->Fill(0.5,(Double_t)fUsePhiWeights||fUsePtWeights||fUseEtaWeights);
2475  // Which event weight was used to weight generating function event-by-event:
2476  if(strcmp(fMultiplicityWeight->Data(),"unit"))
2477  {
2478   fReferenceFlowFlags->Fill(1.5,0.); // 0 = "unit" (default)
2479  } else if(strcmp(fMultiplicityWeight->Data(),"multiplicity"))
2480    {
2481     fReferenceFlowFlags->Fill(1.5,1.); // 1 = "multiplicity"
2482    }
2483
2484 } // end of void AliFlowAnalysisWithCumulants::StoreReferenceFlowFlags()
2485
2486 //================================================================================================================
2487
2488 void AliFlowAnalysisWithCumulants::StoreDiffFlowFlags()
2489 {
2490  // Store all flags for differential flow in profile fDiffFlowFlags.
2491
2492  if(!fDiffFlowFlags)
2493  {
2494   cout<<endl;
2495   cout<<"WARNING: !fDiffFlowFlags is NULL in AFAWC::SRFF() !!!!"<<endl;
2496   cout<<endl;
2497   exit(0);
2498  }
2499
2500  // fDiffFlags->Fill(0.5,(Double_t) ... );
2501
2502 } // end of void AliFlowAnalysisWithCumulants::StoreDiffFlowFlags()
2503
2504 //================================================================================================================
2505
2506 void AliFlowAnalysisWithCumulants::BookAndNestAllLists()
2507 {
2508  // Book and nest all list in base list fHistList.
2509
2510  // a) Book and nest lists for reference flow;
2511  // b) Book and nest lists for differential flow;
2512  // c) Book and nest lists for tuning.
2513
2514  // a) Book and nest all lists for reference flow:
2515  fReferenceFlowList = new TList();
2516  fReferenceFlowList->SetName("Reference Flow");
2517  fReferenceFlowList->SetOwner(kTRUE);
2518  fHistList->Add(fReferenceFlowList);
2519  fReferenceFlowProfiles = new TList();
2520  fReferenceFlowProfiles->SetName("Profiles");
2521  fReferenceFlowProfiles->SetOwner(kTRUE);
2522  fReferenceFlowList->Add(fReferenceFlowProfiles);
2523  fReferenceFlowResults = new TList();
2524  fReferenceFlowResults->SetName("Results");
2525  fReferenceFlowResults->SetOwner(kTRUE);
2526  fReferenceFlowList->Add(fReferenceFlowResults);
2527  // b) Book and nest lists for differential flow:
2528  fDiffFlowList = new TList();
2529  fDiffFlowList->SetName("Differential Flow");
2530  fDiffFlowList->SetOwner(kTRUE);
2531  fHistList->Add(fDiffFlowList);
2532  fDiffFlowProfiles = new TList();
2533  fDiffFlowProfiles->SetName("Profiles");
2534  fDiffFlowProfiles->SetOwner(kTRUE);
2535  fDiffFlowList->Add(fDiffFlowProfiles);
2536  fDiffFlowResults = new TList();
2537  fDiffFlowResults->SetName("Results");
2538  fDiffFlowResults->SetOwner(kTRUE);
2539  fDiffFlowList->Add(fDiffFlowResults);
2540  // c) Book and nest lists for tuning:
2541  if(fTuneParameters)
2542  {
2543   fTuningList = new TList();
2544   fTuningList->SetName("Tuning");
2545   fTuningList->SetOwner(kTRUE);
2546   fHistList->Add(fTuningList);
2547   fTuningProfiles = new TList();
2548   fTuningProfiles->SetName("Profiles");
2549   fTuningProfiles->SetOwner(kTRUE);
2550   fTuningList->Add(fTuningProfiles);
2551   fTuningResults = new TList();
2552   fTuningResults->SetName("Results");
2553   fTuningResults->SetOwner(kTRUE);
2554   fTuningList->Add(fTuningResults);
2555  }
2556
2557 } // end of void AliFlowAnalysisWithCumulants::BookAndNestAllLists()
2558
2559 //================================================================================================================
2560
2561 void AliFlowAnalysisWithCumulants::BookProfileHoldingSettings()
2562 {
2563  // Book profile to hold all analysis settings.
2564
2565  TString analysisSettingsName = "fAnalysisSettings";
2566  fAnalysisSettings = new TProfile(analysisSettingsName.Data(),"Settings for analysis with Generating Function Cumulants",6,0.,6.);
2567  fAnalysisSettings->GetXaxis()->SetLabelSize(0.035);
2568  fAnalysisSettings->GetXaxis()->SetBinLabel(1,"Harmonic");
2569  fAnalysisSettings->Fill(0.5,fHarmonic);
2570  fAnalysisSettings->GetXaxis()->SetBinLabel(2,"Multiple");
2571  fAnalysisSettings->Fill(1.5,fMultiple);
2572  fAnalysisSettings->GetXaxis()->SetBinLabel(3,"r_{0}");
2573  fAnalysisSettings->Fill(2.5,fR0);
2574  fAnalysisSettings->GetXaxis()->SetBinLabel(4,"Print RF results");
2575  fAnalysisSettings->Fill(3.5,fPrintFinalResults[0]);
2576  fAnalysisSettings->GetXaxis()->SetBinLabel(5,"Print RP results");
2577  fAnalysisSettings->Fill(4.5,fPrintFinalResults[1]);
2578  fAnalysisSettings->GetXaxis()->SetBinLabel(6,"Print POI results");
2579  fAnalysisSettings->Fill(5.5,fPrintFinalResults[2]);
2580  fHistList->Add(fAnalysisSettings);
2581
2582 } // end of void AliFlowAnalysisWithCumulants::BookProfileHoldingSettings()
2583
2584 //================================================================================================================
2585
2586 void AliFlowAnalysisWithCumulants::BookCommonHistograms()
2587 {
2588  // Book common control histograms and common histograms for final results.
2589
2590  // Common control histogram:
2591  TString commonHistsName = "AliFlowCommonHistGFC";
2592  fCommonHists = new AliFlowCommonHist(commonHistsName.Data());
2593  fHistList->Add(fCommonHists);
2594  // Common histograms for final results from 2nd order GFC:
2595  TString commonHistResults2ndOrderName = "AliFlowCommonHistResults2ndOrderGFC";
2596  fCommonHistsResults2nd = new AliFlowCommonHistResults(commonHistResults2ndOrderName.Data());
2597  fHistList->Add(fCommonHistsResults2nd);
2598  // Common histograms for final results from 4th order GFC:
2599  TString commonHistResults4thOrderName = "AliFlowCommonHistResults4thOrderGFC";
2600  fCommonHistsResults4th = new AliFlowCommonHistResults(commonHistResults4thOrderName.Data());
2601  fHistList->Add(fCommonHistsResults4th);
2602  // Common histograms for final results from 6th order GFC:
2603  TString commonHistResults6thOrderName = "AliFlowCommonHistResults6thOrderGFC";
2604  fCommonHistsResults6th = new AliFlowCommonHistResults(commonHistResults6thOrderName.Data());
2605  fHistList->Add(fCommonHistsResults6th);
2606  // Common histograms for final results from 8th order GFC:
2607  TString commonHistResults8thOrderName = "AliFlowCommonHistResults8thOrderGFC";
2608  fCommonHistsResults8th = new AliFlowCommonHistResults(commonHistResults8thOrderName.Data());
2609  fHistList->Add(fCommonHistsResults8th);
2610
2611 } // end of void AliFlowAnalysisWithCumulants::BookCommonHistograms()
2612
2613 //================================================================================================================
2614
2615 void AliFlowAnalysisWithCumulants::CheckPointersUsedInMake()
2616 {
2617  // Check pointers used in method Make().
2618
2619  if(!fCommonHists)
2620  {
2621   cout<<endl;
2622   cout<<" WARNING (GFC): !fCommonHists is NULL in CPUIM() !!!!"<<endl;
2623   cout<<endl;
2624   exit(0);
2625  }
2626  if(fUsePhiWeights && !fPhiWeights)
2627  {
2628   cout<<endl;
2629   cout<<" WARNING (GFC): !fPhiWeights is NULL in CPUIM() !!!!"<<endl;
2630   cout<<endl;
2631   exit(0);
2632  }
2633  if(fUsePtWeights && !fPtWeights)
2634  {
2635   cout<<endl;
2636   cout<<" WARNING (GFC): !fPtWeights is NULL in CPUIM() !!!!"<<endl;
2637   cout<<endl;
2638   exit(0);
2639  }
2640  if(fUseEtaWeights && !fEtaWeights)
2641  {
2642   cout<<endl;
2643   cout<<" WARNING (GFC): !fEtaWeights is NULL in CPUIM() !!!!"<<endl;
2644   cout<<endl;
2645   exit(0);
2646  }
2647  if(!fAverageOfSquaredWeight)
2648  {
2649   cout<<endl;
2650   cout<<" WARNING (GFC): !fAverageOfSquaredWeight is NULL in CPUIM() !!!!"<<endl;
2651   cout<<endl;
2652   exit(0);
2653  }
2654  if(!fReferenceFlowGenFun)
2655  {
2656   cout<<endl;
2657   cout<<" WARNING (GFC): !fReferenceFlowGenFun is NULL in CPUIM() !!!!"<<endl;
2658   cout<<endl;
2659   exit(0);
2660  }
2661  if(!fQvectorComponents)
2662  {
2663   cout<<endl;
2664   cout<<" WARNING (GFC): !fQvectorComponents is NULL in CPUIM() !!!!"<<endl;
2665   cout<<endl;
2666   exit(0);
2667  }
2668  if(!fGEBE)
2669  {
2670   cout<<endl;
2671   cout<<"WARNING (GFC): !fGEBE is NULL in CPUIM() !!!!"<<endl;
2672   cout<<endl;
2673   exit(0);
2674  }
2675
2676 } // end of void AliFlowAnalysisWithCumulants::CheckPointersUsedInMake()
2677
2678 //================================================================================================================
2679
2680 void AliFlowAnalysisWithCumulants::CheckPointersUsedInFinish()
2681 {
2682  // Check pointers used in method Finish().
2683
2684  if(!fAnalysisSettings)
2685  {
2686   cout<<endl;
2687   cout<<" WARNING (GFC): fAnalysisSettings is NULL in CPUIF() !!!!"<<endl;
2688   cout<<endl;
2689   exit(0);
2690  }
2691  if(!(fCommonHists && fCommonHists->GetHistMultRP()))
2692  {
2693   cout<<endl;
2694   cout<<" WARNING (GFC): (fCommonHists && fCommonHists->GetHistMultRP) is NULL in CPUIF() !!!!"<<endl;
2695   cout<<endl;
2696   exit(0);
2697  }
2698  if(!fReferenceFlowGenFun)
2699  {
2700   cout<<endl;
2701   cout<<" WARNING (GFC): fReferenceFlowGenFun is NULL in CPUIF() !!!!"<<endl;
2702   cout<<endl;
2703   exit(0);
2704  }
2705  if(!fReferenceFlowCumulants)
2706  {
2707   cout<<endl;
2708   cout<<" WARNING (GFC): fReferenceFlowCumulants is NULL in CPUIF() !!!!"<<endl;
2709   cout<<endl;
2710   exit(0);
2711  }
2712  if(!fQvectorComponents)
2713  {
2714   cout<<endl;
2715   cout<<" WARNING (GFC): fQvectorComponents is NULL in CPUIF() !!!!"<<endl;
2716   cout<<endl;
2717   exit(0);
2718  }
2719  if(!fAverageOfSquaredWeight)
2720  {
2721   cout<<endl;
2722   cout<<" WARNING (GFC): fAverageOfSquaredWeight is NULL in CPUIF() !!!!"<<endl;
2723   cout<<endl;
2724   exit(0);
2725  }
2726  if(!(fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && fCommonHistsResults8th))
2727  {
2728   cout<<endl;
2729   cout<<" WARNING (GFC): fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && "<<endl;
2730   cout<<"                fCommonHistsResults8th is NULL in CPUIF() !!!!"<<endl;
2731   cout<<endl;
2732   exit(0);
2733  }
2734  if(!fReferenceFlow)
2735  {
2736   cout<<endl;
2737   cout<<" WARNING (GFC): fReferenceFlow is NULL in CPUIF() !!!!"<<endl;
2738   cout<<endl;
2739   exit(0);
2740  }
2741  if(!fChi)
2742  {
2743   cout<<endl;
2744   cout<<" WARNING (GFC): fChi is NULL in CPUIF() !!!!"<<endl;
2745   cout<<endl;
2746   exit(0);
2747  }
2748  for(Int_t ri=0;ri<2;ri++)
2749  {
2750   for(Int_t rp=0;rp<2;rp++)
2751   {
2752    for(Int_t pe=0;pe<2;pe++)
2753    {
2754     if(!fDiffFlowGenFun[ri][rp][pe])
2755     {
2756      cout<<endl;
2757      cout<<" WARNING (GFC): "<<Form("fDiffFlowGenFun[%d][%d][%d]",ri,rp,pe)<<" is NULL in CPUIF() !!!!"<<endl;
2758      cout<<endl;
2759      exit(0);
2760     }
2761    }
2762   }
2763  }
2764  for(Int_t rp=0;rp<2;rp++)
2765  {
2766   for(Int_t pe=0;pe<2;pe++)
2767   {
2768    for(Int_t co=0;co<4;co++)
2769    {
2770     if(!fDiffFlowCumulants[rp][pe][co])
2771     {
2772      cout<<endl;
2773      cout<<" WARNING (GFC): "<<Form("fDiffFlowCumulants[%d][%d][%d]",rp,pe,co)<<" is NULL in CPUIF() !!!!"<<endl;
2774      cout<<endl;
2775      exit(0);
2776     }
2777     if(!fDiffFlow[rp][pe][co])
2778     {
2779      cout<<endl;
2780      cout<<" WARNING (GFC): "<<Form("fDiffFlow[%d][%d][%d]",rp,pe,co)<<" is NULL in CPUIF() !!!!"<<endl;
2781      cout<<endl;
2782      exit(0);
2783     }
2784    }
2785   }
2786  }
2787  for(Int_t rp=0;rp<2;rp++)
2788  {
2789   for(Int_t pe=0;pe<2;pe++)
2790   {
2791    if(!fNoOfParticlesInBin[rp][pe])
2792    {
2793     cout<<endl;
2794     cout<<" WARNING (GFC): "<<Form("fNoOfParticlesInBin[%d][%d]",rp,pe)<<" is NULL in CPUIF() !!!!"<<endl;
2795     cout<<endl;
2796     exit(0);
2797    }
2798   }
2799  }
2800
2801 } // end of void AliFlowAnalysisWithCumulants::CheckPointersUsedInFinish()
2802
2803 //================================================================================================================
2804
2805 void AliFlowAnalysisWithCumulants::AccessSettings()
2806 {
2807  // Access the settings for analysis with Generating Function Cumulants.
2808
2809  fHarmonic = (Int_t)fAnalysisSettings->GetBinContent(1);
2810  fMultiple = (Int_t)fAnalysisSettings->GetBinContent(2);
2811  fR0 = (Double_t)fAnalysisSettings->GetBinContent(3);
2812  fPrintFinalResults[0] = (Bool_t)fAnalysisSettings->GetBinContent(4);
2813  fPrintFinalResults[1] = (Bool_t)fAnalysisSettings->GetBinContent(5);
2814  fPrintFinalResults[2] = (Bool_t)fAnalysisSettings->GetBinContent(6);
2815
2816 } // end of AliFlowAnalysisWithCumulants::AccessSettings()
2817
2818 //================================================================================================================
2819
2820 void AliFlowAnalysisWithCumulants::WriteHistograms(TString* outputFileName)
2821 {
2822  // Store the final results in output .root file.
2823
2824  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
2825  //output->WriteObject(fHistList, "cobjGFC","SingleKey");
2826  fHistList->SetName("cobjGFC");
2827  fHistList->SetOwner(kTRUE);
2828  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
2829  delete output;
2830
2831 } // end of void AliFlowAnalysisWithCumulants::WriteHistograms(TString* outputFileName)
2832
2833 //================================================================================================================
2834
2835 void AliFlowAnalysisWithCumulants::WriteHistograms(TString outputFileName)
2836 {
2837  // Store the final results in output .root file.
2838
2839  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
2840  //output->WriteObject(fHistList, "cobjGFC","SingleKey");
2841  fHistList->SetName("cobjGFC");
2842  fHistList->SetOwner(kTRUE);
2843  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
2844  delete output;
2845
2846 } // end of void AliFlowAnalysisWithCumulants::WriteHistograms(TString outputFileName)
2847
2848 //================================================================================================================
2849
2850 void AliFlowAnalysisWithCumulants::WriteHistograms(TDirectoryFile *outputFileName)
2851 {
2852  // Store the final results in output .root file.
2853
2854  fHistList->SetName("cobjGFC");
2855  fHistList->SetOwner(kTRUE);
2856  outputFileName->Add(fHistList);
2857  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
2858
2859 } // end of void AliFlowAnalysisWithCumulants::WriteHistograms(TDirectoryFile *outputFileName)
2860
2861 //================================================================================================================
2862
2863
2864
2865