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 /********************************
  17  * flow analysis with cumulants *
  18  *                              *
  19  * author: Ante Bilandzic       *
  20  *          (anteb@nikhef.nl)   *
  21  *******************************/
  22
  23 #define AliFlowAnalysisWithCumulants_cxx
  24
  25 #include "Riostream.h"
  26 #include "AliFlowCommonConstants.h"
  27 #include "AliFlowCommonHist.h"
  28 #include "AliFlowCommonHistResults.h"
  29 #include "TChain.h"
  30 #include "TFile.h"
  31 #include "TList.h" //NEW
  32 #include "TParticle.h"
  33 #include "TRandom3.h"
  34 #include "TProfile.h"
  35 #include "TProfile2D.h"
  36 #include "TProfile3D.h"
  37 #include "AliFlowEventSimple.h"
  38 #include "AliFlowTrackSimple.h"
  39 #include "AliFlowAnalysisWithCumulants.h"
  40 #include "AliFlowCumuConstants.h"
  41 #include "AliCumulantsFunctions.h"
  42
  43 class TH1;
  44 class TGraph;
  45 class TPave;
  46 class TLatex;
  47 class TMarker;
  48 class TRandom3;
  49 class TObjArray;
  50 class TList;
  51 class TCanvas;
  52 class TSystem;
  53 class TROOT;
  54 class AliFlowVector;
  55 class TVector;
  56
  57 //================================================================================================================
  58
  59 ClassImp(AliFlowAnalysisWithCumulants)
  60
  61 AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants():
  62  fTrack(NULL),
  63  fHistList(NULL),
  64  fWeightsList(NULL),
  65  fR0(0),
  66  fPtMax(0),
  67  fPtMin(0),
  68  fBinWidthPt(0),
  69  fgknBinsPt(0),
  70  fEtaMax(0),
  71  fEtaMin(0),
  72  fBinWidthEta(0),
  73  fgknBinsEta(0),
  74  fAvQx(0),
  75  fAvQy(0),
  76  fAvQ2x(0),
  77  fAvQ2y(0),
  78  fAvMultIntFlowGFC(NULL),
  79  fQVectorComponentsGFC(NULL),
  80  fIntFlowResultsGFC(NULL),
  81  fDiffFlowResults2ndOrderGFC(NULL),
  82  fDiffFlowResults4thOrderGFC(NULL),
  83  fDiffFlowResults6thOrderGFC(NULL),
  84  fDiffFlowResults8thOrderGFC(NULL),
  85  fCommonHistsResults2nd(NULL),
  86  fCommonHistsResults4th(NULL),
  87  fCommonHistsResults6th(NULL),
  88  fCommonHistsResults8th(NULL),
  89  fIntFlowGenFun(NULL),
  90  fIntFlowGenFun4(NULL),//(only for other system of Eq.)
  91  fIntFlowGenFun6(NULL),//(only for other system of Eq.)
  92  fIntFlowGenFun8(NULL),//(only for other system of Eq.)
  93  fIntFlowGenFun16(NULL),//(only for other system of Eq.)
  94  fAvMultIntFlow4GFC(NULL),//(only for other system of Eq.)
  95  fAvMultIntFlow6GFC(NULL),//(only for other system of Eq.)
  96  fAvMultIntFlow8GFC(NULL),//(only for other system of Eq.)
  97  fAvMultIntFlow16GFC(NULL),//(only for other system of Eq.)
  98  fDiffFlowPtRPGenFunRe(NULL),
  99  fDiffFlowPtRPGenFunIm(NULL),
 100  fPtBinRPNoOfParticles(NULL),
 101  fDiffFlowEtaRPGenFunRe(NULL),
 102  fDiffFlowEtaRPGenFunIm(NULL),
 103  fEtaBinRPNoOfParticles(NULL),
 104  fDiffFlowPtPOIGenFunRe(NULL),
 105  fDiffFlowPtPOIGenFunIm(NULL),
 106  fPtBinPOINoOfParticles(NULL),
 107  fDiffFlowEtaPOIGenFunRe(NULL),
 108  fDiffFlowEtaPOIGenFunIm(NULL),
 109  fEtaBinPOINoOfParticles(NULL),
 110  /*
 111  fDiffFlowGenFunRe0(NULL),
 112  fDiffFlowGenFunRe1(NULL),
 113  fDiffFlowGenFunRe2(NULL),
 114  fDiffFlowGenFunRe3(NULL),
 115  fDiffFlowGenFunRe4(NULL),
 116  fDiffFlowGenFunRe5(NULL),
 117  fDiffFlowGenFunRe6(NULL),
 118  fDiffFlowGenFunRe7(NULL),
 119  fDiffFlowGenFunIm0(NULL),
 120  fDiffFlowGenFunIm1(NULL),
 121  fDiffFlowGenFunIm2(NULL),
 122  fDiffFlowGenFunIm3(NULL),
 123  fDiffFlowGenFunIm4(NULL),
 124  fDiffFlowGenFunIm5(NULL),
 125  fDiffFlowGenFunIm6(NULL),
 126  fDiffFlowGenFunIm7(NULL),
 127  */
 128  fCommonHists(NULL),
 129  fOtherEquations(kFALSE),
 130  fUsePhiWeights(kFALSE),
 131  fUsePtWeights(kFALSE),
 132  fUseEtaWeights(kFALSE),
 133  fAverageOfSquaredWeight(NULL)
 134 {
 135  //constructor
 136  fHistList = new TList();
 137  fWeightsList = new TList();
 138  fWeightsList->SetName("Weights");
 139  fR0=AliFlowCumuConstants::fgR0;
 140  //Pt:
 141  fPtMax=AliFlowCommonConstants::GetPtMax();
 142  fPtMin=AliFlowCommonConstants::GetPtMin();
 143  fgknBinsPt=AliFlowCommonConstants::GetNbinsPt();
 144  if(fgknBinsPt)
 145  {
 146   fBinWidthPt=(fPtMax-fPtMin)/fgknBinsPt;
 147  }
 148  //Eta:
 149  fEtaMax=AliFlowCommonConstants::GetEtaMax();
 150  fEtaMin=AliFlowCommonConstants::GetEtaMin();
 151  fgknBinsEta=AliFlowCommonConstants::GetNbinsEta();
 152  if(fgknBinsEta)
 153  {
 154   fBinWidthEta=(fEtaMax-fEtaMin)/fgknBinsEta;
 155  }
 156
 157  fOtherEquations=AliFlowCumuConstants::fgOtherEquations;
 158 }
 159
 160 AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants()
 161 {
 162  //desctructor
 163  delete fHistList;
 164  delete fWeightsList;
 165 }
 166
 167 //================================================================================================================
 168
 169 void AliFlowAnalysisWithCumulants::Init()
 170 {
 171  //various output histograms
 172
 173  //average multiplicity
 174  fAvMultIntFlowGFC = new TProfile("fAvMultIntFlowGFC","Average Weighted Multiplicity",1,0,1,"s");
 175  fAvMultIntFlowGFC->SetXTitle("");
 176  fAvMultIntFlowGFC->SetYTitle("");
 177  fAvMultIntFlowGFC->SetLabelSize(0.06);
 178  fAvMultIntFlowGFC->SetMarkerStyle(25);
 179  fAvMultIntFlowGFC->SetLabelOffset(0.01);
 180  (fAvMultIntFlowGFC->GetXaxis())->SetBinLabel(1,"Average Weighted Multiplicity");
 181  fHistList->Add(fAvMultIntFlowGFC);
 182
 183  //averages of Q-vector components (1st bin: <Q_x>, 2nd bin: <Q_y>, 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>)
 184  fQVectorComponentsGFC = new TProfile("fQVectorComponentsGFC","Average of Q-vector components",4,0.,4.);
 185  fQVectorComponentsGFC->SetXTitle("");
 186  fQVectorComponentsGFC->SetYTitle("");
 187  fQVectorComponentsGFC->SetLabelSize(0.06);
 188  fQVectorComponentsGFC->SetMarkerStyle(25);
 189  (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(1,"<Q_{x}>");
 190  (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(2,"<Q_{y}>");
 191  (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(3,"<Q_{x}^{2}>");
 192  (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(4,"<Q_{y}^{2}>");
 193  fHistList->Add(fQVectorComponentsGFC);
 194
 195  //final results for integrated flow (v_n{2}, v_n{4},..., v_n{16}) from cumulants (by default n=2)
 196  fIntFlowResultsGFC = new TH1D("fIntFlowResultsGFC","Integrated Flow From Cumulants (Generating Function)",4,0,4);
 197  fIntFlowResultsGFC->SetXTitle("");
 198  fIntFlowResultsGFC->SetYTitle("");
 199  fIntFlowResultsGFC->SetLabelSize(0.06);
 200  //fIntFlowResultsGFC->SetTickLength(1);
 201  fIntFlowResultsGFC->SetMarkerStyle(25);
 202  (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(1,"v_{n}{2}");
 203  (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(2,"v_{n}{4}");
 204  (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(3,"v_{n}{6}");
 205  (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(4,"v_{n}{8}");
 206  fHistList->Add(fIntFlowResultsGFC);
 207
 208  //final results for differential flow v_p/n{2} (by default p=n=2)
 209  fDiffFlowResults2ndOrderGFC = new TH1D("fDiffFlowResults2ndOrderGFC","v'_2/2{2}",fgknBinsPt,fPtMin,fPtMax);
 210  fDiffFlowResults2ndOrderGFC->SetXTitle("pt [GeV]");
 211  fDiffFlowResults2ndOrderGFC->SetYTitle("");
 212  fHistList->Add(fDiffFlowResults2ndOrderGFC);
 213
 214  //final results for differential flow v_p/n{4} (by default p=n=2)
 215  fDiffFlowResults4thOrderGFC = new TH1D("fDiffFlowResults4thOrderGFC","v'_2/2{4}",fgknBinsPt,fPtMin,fPtMax);
 216  fDiffFlowResults4thOrderGFC->SetXTitle("pt [GeV]");
 217  fDiffFlowResults4thOrderGFC->SetYTitle("");
 218  fHistList->Add(fDiffFlowResults4thOrderGFC);
 219
 220  //final results for differential flow v_p/n{6} (by default p=n=2)
 221  fDiffFlowResults6thOrderGFC = new TH1D("fDiffFlowResults6thOrderGFC","v'_2/2{6}",fgknBinsPt,fPtMin,fPtMax);
 222  fDiffFlowResults6thOrderGFC->SetXTitle("pt [GeV]");
 223  fDiffFlowResults6thOrderGFC->SetYTitle("");
 224  fHistList->Add(fDiffFlowResults6thOrderGFC);
 225
 226  //final results for differential flow v_p/n{8} (by default p=n=2)
 227  fDiffFlowResults8thOrderGFC = new TH1D("fDiffFlowResults8thOrderGFC","v'_2/2{8}",fgknBinsPt,fPtMin,fPtMax);
 228  fDiffFlowResults8thOrderGFC->SetXTitle("pt [GeV]");
 229  fDiffFlowResults8thOrderGFC->SetYTitle("");
 230  fHistList->Add(fDiffFlowResults8thOrderGFC);
 231
 232  //avarage of the generating function for integrated flow <G[p][q]>
 233  fIntFlowGenFun = new TProfile2D("fIntFlowGenFun","<G[p][q]>",fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 234  fIntFlowGenFun->SetXTitle("p");
 235  fIntFlowGenFun->SetYTitle("q");
 236  fHistList->Add(fIntFlowGenFun);
 237
 238  if(fOtherEquations)
 239  {
 240   //avarage of the generating function for integrated flow <G[p][q]> (only for other system of Eq. - up to 4th order)
 241   fIntFlowGenFun4 = new TProfile2D("fIntFlowGenFun4","<G4[p4][q4]>",fgkPmax4,0.,(Double_t)fgkPmax4,fgkQmax4,0.,(Double_t)fgkQmax4);
 242   fIntFlowGenFun4->SetXTitle("p4");
 243   fIntFlowGenFun4->SetYTitle("q4");
 244   fHistList->Add(fIntFlowGenFun4);
 245
 246   //avarage of the generating function for integrated flow <G[p][q]> (only for other system of Eq. - up to 6th order)
 247   fIntFlowGenFun6 = new TProfile2D("fIntFlowGenFun6","<G6[p6][q6]>",fgkPmax6,0.,(Double_t)fgkPmax6,fgkQmax6,0.,(Double_t)fgkQmax6);
 248   fIntFlowGenFun6->SetXTitle("p6");
 249   fIntFlowGenFun6->SetYTitle("q6");
 250   fHistList->Add(fIntFlowGenFun6);
 251
 252   //avarage of the generating function for integrated flow <G[p][q]> (only for other system of Eq. - up to 8th order)
 253   fIntFlowGenFun8 = new TProfile2D("fIntFlowGenFun8","<G8[p8][q8]>",fgkPmax8,0.,(Double_t)fgkPmax8,fgkQmax8,0.,(Double_t)fgkQmax8);
 254   fIntFlowGenFun8->SetXTitle("p8");
 255   fIntFlowGenFun8->SetYTitle("q8");
 256   fHistList->Add(fIntFlowGenFun8);
 257
 258   //avarage of the generating function for integrated flow <G[p][q]> (only for other system of Eq. - up to 16th order)
 259   fIntFlowGenFun16 = new TProfile2D("fIntFlowGenFun16","<G16[p16][q16]>",fgkPmax16,0.,(Double_t)fgkPmax16,fgkQmax16,0.,(Double_t)fgkQmax16);
 260   fIntFlowGenFun16->SetXTitle("p16");
 261   fIntFlowGenFun16->SetYTitle("q16");
 262   fHistList->Add(fIntFlowGenFun16);
 263
 264   //average multiplicity (only for other system of Eq. - up to 4th order)
 265   fAvMultIntFlow4GFC = new TProfile("fAvMultIntFlow4GFC","Average Multiplicity",1,0,1,"s");
 266   fAvMultIntFlow4GFC->SetXTitle("");
 267   fAvMultIntFlow4GFC->SetYTitle("");
 268   fAvMultIntFlow4GFC->SetLabelSize(0.06);
 269   fAvMultIntFlow4GFC->SetMarkerStyle(25);
 270   fAvMultIntFlow4GFC->SetLabelOffset(0.01);
 271   (fAvMultIntFlow4GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity");
 272   fHistList->Add(fAvMultIntFlow4GFC);
 273
 274   //average multiplicity (only for other system of Eq. - up to 6th order)
 275   fAvMultIntFlow6GFC = new TProfile("fAvMultIntFlow6GFC","Average Multiplicity",1,0,1,"s");
 276   fAvMultIntFlow6GFC->SetXTitle("");
 277   fAvMultIntFlow6GFC->SetYTitle("");
 278   fAvMultIntFlow6GFC->SetLabelSize(0.06);
 279   fAvMultIntFlow6GFC->SetMarkerStyle(25);
 280   fAvMultIntFlow6GFC->SetLabelOffset(0.01);
 281   (fAvMultIntFlow6GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity");
 282   fHistList->Add(fAvMultIntFlow6GFC);
 283
 284   //average multiplicity (only for other system of Eq. - up to 8th order)
 285   fAvMultIntFlow8GFC = new TProfile("fAvMultIntFlow8GFC","Average Multiplicity",1,0,1,"s");
 286   fAvMultIntFlow8GFC->SetXTitle("");
 287   fAvMultIntFlow8GFC->SetYTitle("");
 288   fAvMultIntFlow8GFC->SetLabelSize(0.06);
 289   fAvMultIntFlow8GFC->SetMarkerStyle(25);
 290   fAvMultIntFlow8GFC->SetLabelOffset(0.01);
 291   (fAvMultIntFlow8GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity");
 292   fHistList->Add(fAvMultIntFlow8GFC);
 293
 294   //average multiplicity (only for other system of Eq. - up to 16th order)
 295   fAvMultIntFlow16GFC = new TProfile("fAvMultIntFlow16GFC","Average Multiplicity",1,0,1,"s");
 296   fAvMultIntFlow16GFC->SetXTitle("");
 297   fAvMultIntFlow16GFC->SetYTitle("");
 298   fAvMultIntFlow16GFC->SetLabelSize(0.06);
 299   fAvMultIntFlow16GFC->SetMarkerStyle(25);
 300   fAvMultIntFlow16GFC->SetLabelOffset(0.01);
 301   (fAvMultIntFlow16GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity");
 302   fHistList->Add(fAvMultIntFlow16GFC);
 303  }
 304
 305  //avarage of the real part of generating function for differential flow in Pt <Re(D[b][p][q])>
 306  fDiffFlowPtRPGenFunRe = new TProfile3D("fDiffFlowPtRPGenFunRe","<Re(D[b][p][q])>",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 307  fDiffFlowPtRPGenFunRe->SetXTitle("b");
 308  fDiffFlowPtRPGenFunRe->SetYTitle("p");
 309  fDiffFlowPtRPGenFunRe->SetZTitle("q");
 310  fDiffFlowPtRPGenFunRe->SetTitleOffset(1.44,"X");
 311  fDiffFlowPtRPGenFunRe->SetTitleOffset(1.44,"Y");
 312  fHistList->Add(fDiffFlowPtRPGenFunRe);
 313
 314  //avarage of the imaginary part of generating function for differential flow in Pt <Im(D[b][p][q])>
 315  fDiffFlowPtRPGenFunIm = new TProfile3D("fDiffFlowPtRPGenFunIm","<Im(D[b][p][q])>",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 316  fDiffFlowPtRPGenFunIm->SetXTitle("b");
 317  fDiffFlowPtRPGenFunIm->SetYTitle("p");
 318  fDiffFlowPtRPGenFunIm->SetZTitle("q");
 319  fDiffFlowPtRPGenFunIm->SetTitleOffset(1.44,"X");
 320  fDiffFlowPtRPGenFunIm->SetTitleOffset(1.44,"Y");
 321  fHistList->Add(fDiffFlowPtRPGenFunIm);
 322
 323  //number of particles per pt bin
 324  fPtBinRPNoOfParticles = new TProfile("fPtBinRPNoOfParticles","Number of particles per #p_{t} bin",fgknBinsPt,fPtMin,fPtMax);
 325  fPtBinRPNoOfParticles->SetXTitle("pt [GeV]");
 326  fPtBinRPNoOfParticles->SetYTitle("");
 327  fHistList->Add(fPtBinRPNoOfParticles);
 328
 329  //avarage of the real part of generating function for differential flow in Eta <Re(D[b][p][q])>
 330  fDiffFlowEtaRPGenFunRe = new TProfile3D("fDiffFlowEtaRPGenFunRe","<Re(D[b][p][q])>",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 331  fDiffFlowEtaRPGenFunRe->SetXTitle("b");
 332  fDiffFlowEtaRPGenFunRe->SetYTitle("p");
 333  fDiffFlowEtaRPGenFunRe->SetZTitle("q");
 334  fDiffFlowEtaRPGenFunRe->SetTitleOffset(1.44,"X");
 335  fDiffFlowEtaRPGenFunRe->SetTitleOffset(1.44,"Y");
 336  fHistList->Add(fDiffFlowEtaRPGenFunRe);
 337
 338  //avarage of the imaginary part of generating function for differential flow in Eta <Im(D[b][p][q])>
 339  fDiffFlowEtaRPGenFunIm = new TProfile3D("fDiffFlowEtaRPGenFunIm","<Im(D[b][p][q])>",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 340  fDiffFlowEtaRPGenFunIm->SetXTitle("b");
 341  fDiffFlowEtaRPGenFunIm->SetYTitle("p");
 342  fDiffFlowEtaRPGenFunIm->SetZTitle("q");
 343  fDiffFlowEtaRPGenFunIm->SetTitleOffset(1.44,"X");
 344  fDiffFlowEtaRPGenFunIm->SetTitleOffset(1.44,"Y");
 345  fHistList->Add(fDiffFlowEtaRPGenFunIm);
 346
 347  //number of particles per eta bin
 348  fEtaBinRPNoOfParticles = new TProfile("fEtaBinRPNoOfParticles","Number of particles per #eta bin",fgknBinsEta,fEtaMin,fEtaMax);
 349  fEtaBinRPNoOfParticles->SetXTitle("#eta");
 350  fEtaBinRPNoOfParticles->SetYTitle("");
 351  fHistList->Add(fEtaBinRPNoOfParticles);
 352
 353  //avarage of the real part of generating function for differential flow in Pt <Re(D[b][p][q])>
 354  fDiffFlowPtPOIGenFunRe = new TProfile3D("fDiffFlowPtPOIGenFunRe","<Re(D[b][p][q])>",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 355  fDiffFlowPtPOIGenFunRe->SetXTitle("b");
 356  fDiffFlowPtPOIGenFunRe->SetYTitle("p");
 357  fDiffFlowPtPOIGenFunRe->SetZTitle("q");
 358  fDiffFlowPtPOIGenFunRe->SetTitleOffset(1.44,"X");
 359  fDiffFlowPtPOIGenFunRe->SetTitleOffset(1.44,"Y");
 360  fHistList->Add(fDiffFlowPtPOIGenFunRe);
 361
 362  //avarage of the imaginary part of generating function for differential flow in Pt <Im(D[b][p][q])>
 363  fDiffFlowPtPOIGenFunIm = new TProfile3D("fDiffFlowPtPOIGenFunIm","<Im(D[b][p][q])>",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 364  fDiffFlowPtPOIGenFunIm->SetXTitle("b");
 365  fDiffFlowPtPOIGenFunIm->SetYTitle("p");
 366  fDiffFlowPtPOIGenFunIm->SetZTitle("q");
 367  fDiffFlowPtPOIGenFunIm->SetTitleOffset(1.44,"X");
 368  fDiffFlowPtPOIGenFunIm->SetTitleOffset(1.44,"Y");
 369  fHistList->Add(fDiffFlowPtPOIGenFunIm);
 370
 371  //number of particles per pt bin
 372  fPtBinPOINoOfParticles = new TProfile("fPtBinPOINoOfParticles","Number of particles per #p_{t} bin",fgknBinsPt,fPtMin,fPtMax);
 373  fPtBinPOINoOfParticles->SetXTitle("pt [GeV]");
 374  fPtBinPOINoOfParticles->SetYTitle("");
 375  fHistList->Add(fPtBinPOINoOfParticles);
 376
 377  //avarage of the real part of generating function for differential flow in Eta <Re(D[b][p][q])>
 378  fDiffFlowEtaPOIGenFunRe = new TProfile3D("fDiffFlowEtaPOIGenFunRe","<Re(D[b][p][q])>",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 379  fDiffFlowEtaPOIGenFunRe->SetXTitle("b");
 380  fDiffFlowEtaPOIGenFunRe->SetYTitle("p");
 381  fDiffFlowEtaPOIGenFunRe->SetZTitle("q");
 382  fDiffFlowEtaPOIGenFunRe->SetTitleOffset(1.44,"X");
 383  fDiffFlowEtaPOIGenFunRe->SetTitleOffset(1.44,"Y");
 384  fHistList->Add(fDiffFlowEtaPOIGenFunRe);
 385
 386  //avarage of the imaginary part of generating function for differential flow in Eta <Im(D[b][p][q])>
 387  fDiffFlowEtaPOIGenFunIm = new TProfile3D("fDiffFlowEtaPOIGenFunIm","<Im(D[b][p][q])>",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax);
 388  fDiffFlowEtaPOIGenFunIm->SetXTitle("b");
 389  fDiffFlowEtaPOIGenFunIm->SetYTitle("p");
 390  fDiffFlowEtaPOIGenFunIm->SetZTitle("q");
 391  fDiffFlowEtaPOIGenFunIm->SetTitleOffset(1.44,"X");
 392  fDiffFlowEtaPOIGenFunIm->SetTitleOffset(1.44,"Y");
 393  fHistList->Add(fDiffFlowEtaPOIGenFunIm);
 394
 395  //number of particles per eta bin
 396  fEtaBinPOINoOfParticles = new TProfile("fEtaBinPOINoOfParticles","Number of particles per #eta bin",fgknBinsEta,fEtaMin,fEtaMax);
 397  fEtaBinPOINoOfParticles->SetXTitle("#eta");
 398  fEtaBinPOINoOfParticles->SetYTitle("");
 399  fHistList->Add(fEtaBinPOINoOfParticles);
 400
 401  /*
 402  fDiffFlowGenFunRe0 = new TProfile2D("fDiffFlowGenFunRe0","Re(<D[b][0][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 403  fDiffFlowGenFunRe0->SetXTitle("b");
 404  fDiffFlowGenFunRe0->SetYTitle("q");
 405  fHistList->Add(fDiffFlowGenFunRe0);
 406
 407  fDiffFlowGenFunIm0 = new TProfile2D("fDiffFlowGcout<<"HEY M1"<<endl;enFunIm0","Im(<D[b][0][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 408  fDiffFlowGenFunIm0->SetXTitle("b");
 409  fDiffFlowGenFunIm0->SetYTitle("q");
 410  fHistList->Add(fDiffFlowGenFunIm0);
 411
 412  fDiffFlowGenFunRe1 = new TProfile2D("fDiffFlowGenFunRe1","Re(<D[b][1][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 413  fDiffFlowGenFunRe1->SetXTitle("b");
 414  fDiffFlowGenFunRe1->SetYTitle("q");
 415  fHistList->Add(fDiffFlowGenFunRe1);
 416
 417  fDiffFlowGenFunIm1 = new TProfile2D("fDiffFlowGenFunIm1","Im(<D[b][1][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 418  fDiffFlowGenFunIm1->SetXTitle("b");
 419  fDiffFlowGenFunIm1->SetYTitle("q");
 420  fHistList->Add(fDiffFlowGenFunIm1);
 421
 422  fDiffFlowGenFunRe2 = new TProfile2D("fDiffFlowGenFunRe2","Re(<D[b][2][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 423  fDiffFlowGenFunRe2->SetXTitle("b");
 424  fDiffFlowGenFunRe2->SetYTitle("q");
 425  fHistList->Add(fDiffFlowGenFunRe2);
 426
 427  fDiffFlowGenFunIm2 = new TProfile2D("fDiffFlowGenFunIm2","Im(<D[b][2][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 428  fDiffFlowGenFunIm2->SetXTitle("b");
 429  fDiffFlowGenFunIm2->SetYTitle("q");
 430  fHistList->Add(fDiffFlowGenFunIm2);
 431
 432  fDiffFlowGenFunRe3 = new TProfile2D("fDiffFlowGenFunRe3","Re(<D[b][3][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 433  fDiffFlowGenFunRe3->SetXTitle("b");
 434  fDiffFlowGenFunRe3->SetYTitle("q");
 435  fHistList->Add(fDiffFlowGenFunRe3);
 436
 437  fDiffFlowGenFunIm3 = new TProfile2D("fDiffFlowGenFunIm3","Im(<D[b][3][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 438  fDiffFlowGenFunIm3->SetXTitle("b");
 439  fDiffFlowGenFunIm3->SetYTitle("q");
 440  fHistList->Add(fDiffFlowGenFunIm3);
 441
 442  fDiffFlowGenFunRe4 = new TProfile2D("fDiffFlowGenFunRe4","Re(<D[b][4][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 443  fDiffFlowGenFunRe4->SetXTitle("b");
 444  fDiffFlowGenFunRe4->SetYTitle("q");
 445  fHistList->Add(fDiffFlowGenFunRe4);
 446
 447  fDiffFlowGenFunIm4 = new TProfile2D("fDiffFlowGenFunIm4","Im(<D[b][4][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 448  fDiffFlowGenFunIm4->SetXTitle("b");
 449  fDiffFlowGenFunIm4->SetYTitle("q");
 450  fHistList->Add(fDiffFlowGenFunIm4);
 451
 452  fDiffFlowGenFunRe5 = new TProfile2D("fDiffFlowGenFunRe5","Re(<D[b][5][q]>)",fgkQmax,0.,(Double_t)fgkQmax,fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth));
 453  fDiffFlowGenFunRe5->SetXTitle("b");
 454  fDiffFlowGenFunRe5->SetYTitle("q");
 455  fHistList->Add(fDiffFlowGenFunRe5);
 456
 457  fDiffFlowGenFunIm5 = new TProfile2D("fDiffFlowGenFunIm5","Im(<D[b][5][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 458  fDiffFlowGenFunIm5->SetXTitle("b");
 459  fDiffFlowGenFunIm5->SetYTitle("q");
 460  fHistList->Add(fDiffFlowGenFunIm5);
 461
 462  fDiffFlowGenFunRe6 = new TProfile2D("fDiffFlowGenFunRe6","Re(<D[b][6][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 463  fDiffFlowGenFunRe6->SetXTitle("b");
 464  fDiffFlowGenFunRe6->SetYTitle("q");
 465  fHistList->Add(fDiffFlowGenFunRe6);
 466
 467  fDiffFlowGenFunIm6 = new TProfile2D("fDiffFlowGenFunIm6","Im(<D[b][6][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 468  fDiffFlowGenFunIm6->SetXTitle("b");
 469  fDiffFlowGenFunIm6->SetYTitle("q");
 470  fHistList->Add(fDiffFlowGenFunIm6);
 471
 472  fDiffFlowGenFunRe7 = new TProfile2D("fDiffFlowGenFunRe7","Re(<D[b][7][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 473  fDiffFlowGenFunRe7->SetXTitle("b");
 474  fDiffFlowGenFunRe7->SetYTitle("q");
 475  fHistList->Add(fDiffFlowGenFunRe7);
 476
 477  fDiffFlowGenFunIm7 = new TProfile2D("fDiffFlowGenFunIm7","Im(<D[b][7][q]>)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax);
 478  fDiffFlowGenFunIm7->SetXTitle("b");
 479  fDiffFlowGenFunIm7->SetYTitle("q");
 480  fHistList->Add(fDiffFlowGenFunIm7);
 481  */
 482
 483  //common control histograms
 484  fCommonHists = new AliFlowCommonHist("AliFlowCommonHistGFC");
 485  fHistList->Add(fCommonHists);
 486
 487  //common histograms for final results (2nd order)
 488  fCommonHistsResults2nd = new AliFlowCommonHistResults("AliFlowCommonHistResults2ndOrderGFC");
 489  fHistList->Add(fCommonHistsResults2nd);
 490
 491  //common histograms for final results (4th order)
 492  fCommonHistsResults4th = new AliFlowCommonHistResults("AliFlowCommonHistResults4thOrderGFC");
 493  fHistList->Add(fCommonHistsResults4th);
 494
 495  //common histograms for final results (6th order)
 496  fCommonHistsResults6th = new AliFlowCommonHistResults("AliFlowCommonHistResults6thOrderGFC");
 497  fHistList->Add(fCommonHistsResults6th);
 498
 499  //common histograms for final results (8th order)
 500  fCommonHistsResults8th = new AliFlowCommonHistResults("AliFlowCommonHistResults8thOrderGFC");
 501  fHistList->Add(fCommonHistsResults8th);
 502
 503  //<w^2>
 504  fAverageOfSquaredWeight = new TProfile("fAverageOfSquaredWeight","<w^{2}>",1,0,1);
 505  fAverageOfSquaredWeight->SetLabelSize(0.06);
 506  fAverageOfSquaredWeight->SetMarkerStyle(25);
 507  fAverageOfSquaredWeight->SetLabelOffset(0.01);
 508  (fAverageOfSquaredWeight->GetXaxis())->SetBinLabel(1,"<w^{2}>");
 509  fHistList->Add(fAverageOfSquaredWeight);
 510
 511  // add list fWeightsList with weights to the main list
 512  fHistList->Add(fWeightsList);
 513 }//end of Init()
 514
 515 //================================================================================================================
 516
 517 void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent)
 518 {
 519  //running over data:
 520  Int_t nPrim = anEvent->NumberOfTracks(); //total multiplicity
 521
 522  Int_t nEventNSelTracksIntFlow = anEvent->GetEventNSelTracksIntFlow(); //selected multiplicity (particles used for int. flow)
 523
 524  Int_t n = 2; // int flow harmonic (to be improved)
 525
 526  //---------------------------------------------------------------------------------------------------------
 527  // weights:
 528  Bool_t useWeights = fUsePhiWeights||fUsePtWeights||fUseEtaWeights;
 529
 530  TH1F *phiWeights = NULL; // histogram with phi weights
 531  TH1D *ptWeights  = NULL; // histogram with pt weights
 532  TH1D *etaWeights = NULL; // histogram with eta weights
 533
 534  Double_t wPhi = 1.; // phi weight
 535  Double_t wPt  = 1.; // pt weight
 536  Double_t wEta = 1.; // eta weight
 537
 538  if(useWeights)
 539  {
 540   if(!fWeightsList)
 541   {
 542    cout<<" WARNING: fWeightsList is NULL pointer. "<<endl;
 543    exit(0);
 544   }
 545   if(fUsePhiWeights)
 546   {
 547    phiWeights = dynamic_cast<TH1F *>(fWeightsList->FindObject("phi_weights"));
 548    if(!phiWeights)
 549    {
 550     cout<<" WARNING: couldn't access the histogram with phi weights. "<<endl;
 551     exit(0);
 552    }
 553   }
 554   if(fUsePtWeights)
 555   {
 556    ptWeights = dynamic_cast<TH1D *>(fWeightsList->FindObject("pt_weights"));
 557    if(!ptWeights)
 558    {
 559     cout<<" WARNING: couldn't access the histogram with pt weights. "<<endl;
 560     exit(0);
 561    }
 562   }
 563   if(fUseEtaWeights)
 564   {
 565    etaWeights = dynamic_cast<TH1D *>(fWeightsList->FindObject("eta_weights"));
 566    if(!etaWeights)
 567    {
 568     cout<<" WARNING: couldn't access the histogram with eta weights. "<<endl;
 569     exit(0);
 570    }
 571   }
 572  }
 573  //---------------------------------------------------------------------------------------------------------
 574
 575  if(nEventNSelTracksIntFlow>9) //generating function formalism applied here make sense only for selected multiplicity >= 10
 576  {
 577  //fill the common control histograms
 578  fCommonHists->FillControlHistograms(anEvent);
 579
 580  //initializing the generating function G[p][q] for integrated flow
 581  Double_t genfunG[fgkPmax][fgkQmax];
 582
 583  for(Int_t p=0;p<fgkPmax;p++)
 584   {
 585    for(Int_t q=0;q<fgkQmax;q++)
 586    {
 587     genfunG[p][q]=1.;
 588    }
 589   }
 590
 591  Int_t nSelTracksIntFlow = 0; //cross-checking the selected multiplicity
 592
 593  Double_t dPhi = 0.;
 594  Double_t dPt  = 0.;
 595  Double_t dEta = 0.;
 596  Int_t nBinsPhi = 0;
 597
 598  for(Int_t i=0;i<nPrim;i++)
 599  {
 600   fTrack=anEvent->GetTrack(i);
 601   if(fTrack && fTrack->UseForIntegratedFlow())
 602   {
 603    nSelTracksIntFlow++;
 604    // get azimuthal angle, momentum and pseudorapidity of a particle:
 605    dPhi = fTrack->Phi();
 606    dPt  = fTrack->Pt();
 607    dEta = fTrack->Eta();
 608    // phi weights:
 609    if(fUsePhiWeights)
 610    {
 611     nBinsPhi = phiWeights->GetNbinsX();
 612     if(nBinsPhi)
 613     {
 614      wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
 615     }
 616    }
 617    // pt weights:
 618    if(fUsePtWeights)
 619    {
 620     if(fBinWidthPt)
 621     {
 622      wPt = ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fBinWidthPt)));
 623     }
 624    }
 625    // eta weights:
 626    if(fUseEtaWeights)
 627    {
 628     if(fBinWidthEta)
 629     {
 630      wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fBinWidthEta)));
 631     }
 632    }
 633    // evaluate the generating function:
 634    for(Int_t p=0;p<fgkPmax;p++)
 635    {
 636     for(Int_t q=0;q<fgkQmax;q++)
 637     {
 638      genfunG[p][q]*=(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nEventNSelTracksIntFlow)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax));
 639     }
 640    }
 641    // calculate <w^2>
 642    fAverageOfSquaredWeight->Fill(0.,pow(wPhi*wPt*wEta,2.),1);
 643   }
 644  } // end of for(Int_t i=0;i<nPrim;i++)
 645
 646  //storing the value of G[p][q] in 2D profile in order to get automatically the avarage <G[p][q]>
 647  for(Int_t p=0;p<fgkPmax;p++)
 648  {
 649   for(Int_t q=0;q<fgkQmax;q++)
 650   {
 651    fIntFlowGenFun->Fill((Double_t)p,(Double_t)q,genfunG[p][q],1);
 652   }
 653  }
 654
 655  //storing the selected multiplicity (if fAvMultIntFlow is not filled here then you had wrongly selected the particles used for integrated flow)
 656  if(nSelTracksIntFlow==nEventNSelTracksIntFlow)
 657  {
 658   fAvMultIntFlowGFC->Fill(0.,nSelTracksIntFlow,1);
 659  }
 660
 661  // calculating Q-vector of event (needed for errors)
 662  AliFlowVector fQVector;
 663  fQVector.Set(0.,0.);
 664  fQVector.SetMult(0);
 665  fQVector=anEvent->GetQ(1*n,fWeightsList,fUsePhiWeights,fUsePtWeights,fUseEtaWeights); // get the Q vector for this event
 666  fQVectorComponentsGFC->Fill(0.,fQVector.X(),1);         // in the 1st bin fill Q_x
 667  fQVectorComponentsGFC->Fill(1.,fQVector.Y(),1);         // in the 2nd bin fill Q_y
 668  fQVectorComponentsGFC->Fill(2.,pow(fQVector.X(),2.),1); // in the 3rd bin fill (Q_x)^2
 669  fQVectorComponentsGFC->Fill(3.,pow(fQVector.Y(),2.),1); // in the 4th bin fill (Q_y)^2
 670
 671  //3D profiles for differential flow in pt and eta
 672  //second loop over data: evaluating the generating function D[b][p][q] for differential flow
 673  //remark 0: D[b][p][q] is a complex number => real and imaginary part are calculated separately
 674  //remark 1: note that bellow G[p][q] is needed, the value of generating function for integrated flow for the CURRENT event
 675  //remark 2: results are stored in 3D profiles in order to automatically get <Re(D[b][p][q])> and <Im(D[b][p][q])>
 676  for(Int_t i=0;i<nPrim;i++)
 677  {
 678   fTrack=anEvent->GetTrack(i);
 679   if(fTrack)
 680   {
 681    if(fTrack->UseForIntegratedFlow() && fTrack->UseForDifferentialFlow())
 682    {
 683     fPtBinPOINoOfParticles->Fill(fTrack->Pt(),fTrack->Pt(),1.);
 684     fEtaBinPOINoOfParticles->Fill(fTrack->Eta(),fTrack->Eta(),1.);
 685     for(Int_t p=0;p<fgkPmax;p++)
 686     {
 687      for(Int_t q=0;q<fgkQmax;q++)
 688      {
 689       //real part (Pt)
 690       fDiffFlowPtPOIGenFunRe->Fill(fTrack->Pt()/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 691       //imaginary part (Pt)
 692       fDiffFlowPtPOIGenFunIm->Fill(fTrack->Pt()/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 693       //real part (Eta)
 694       fDiffFlowEtaPOIGenFunRe->Fill(fTrack->Eta()/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 695       //imaginary part (Eta)
 696       fDiffFlowEtaPOIGenFunIm->Fill(fTrack->Eta()/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 697      }
 698     }
 699    }
 700    else if(fTrack->UseForDifferentialFlow() && !(fTrack->UseForIntegratedFlow()))
 701    {
 702     fPtBinPOINoOfParticles->Fill(fTrack->Pt(),fTrack->Pt(),1.);
 703     fEtaBinPOINoOfParticles->Fill(fTrack->Eta(),fTrack->Eta(),1.);
 704     for(Int_t p=0;p<fgkPmax;p++)
 705     {
 706      for(Int_t q=0;q<fgkQmax;q++)
 707      {
 708       //real part (Pt)
 709       fDiffFlowPtPOIGenFunRe->Fill(fTrack->Pt()/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi()),1.);
 710       //imaginary part (Pt)
 711       fDiffFlowPtPOIGenFunIm->Fill(fTrack->Pt()/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi()),1.);
 712       //real part (Eta)
 713       fDiffFlowEtaPOIGenFunRe->Fill(fTrack->Eta()/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi()),1.);
 714       //imaginary part (Eta)
 715       fDiffFlowEtaPOIGenFunIm->Fill(fTrack->Eta()/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi()),1.);
 716      }
 717     }
 718    }
 719    //RP:
 720    if(fTrack->UseForIntegratedFlow())
 721    {
 722     fPtBinRPNoOfParticles->Fill(fTrack->Pt(),fTrack->Pt(),1.);
 723     fEtaBinRPNoOfParticles->Fill(fTrack->Eta(),fTrack->Eta(),1.);
 724     for(Int_t p=0;p<fgkPmax;p++)
 725     {
 726      for(Int_t q=0;q<fgkQmax;q++)
 727      {
 728       //real part (Pt)
 729       fDiffFlowPtRPGenFunRe->Fill(fTrack->Pt()/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 730       //imaginary part (Pt)
 731       fDiffFlowPtRPGenFunIm->Fill(fTrack->Pt()/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 732       //real part (Eta)
 733       fDiffFlowEtaRPGenFunRe->Fill(fTrack->Eta()/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 734       //imaginary part (Eta)
 735       fDiffFlowEtaRPGenFunIm->Fill(fTrack->Eta()/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(p+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 736      }
 737     }
 738    }//end of if(fTrack->UseForIntegratedFlow())
 739   }//end of if(fTrack)
 740  }//ending the second loop over data
 741
 742
 743
 744  /*
 745  //sixteen 2D profiles for differential flow
 746  for(Int_t i=0;i<nPrim;i++){
 747   fTrack=anEvent->GetTrack(i);
 748   if (fTrack && fTrack->UseForDifferentialFlow()){
 749    //for(Int_t p=0;p<fgkPmax;p++){
 750     for(Int_t q=0;q<fgkQmax;q++){
 751      //real part
 752      fDiffFlowGenFunRe0->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[0][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(0.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 753      //imaginary part
 754      fDiffFlowGenFunIm0->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[0][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(0.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 755      //-----------------------------------------------------------------------
 756      //real part
 757      fDiffFlowGenFunRe1->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[1][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(1.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 758      //imaginary part
 759      fDiffFlowGenFunIm1->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[1][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(1.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 760     //-----------------------------------------------------------------------
 761      //real part
 762      fDiffFlowGenFunRe2->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[2][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(2.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 763      //imaginary part
 764      fDiffFlowGenFunIm2->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[2][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(2.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 765      //-----------------------------------------------------------------------
 766      //real part
 767      fDiffFlowGenFunRe3->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[3][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(3.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 768      //imaginary part
 769      fDiffFlowGenFunIm3->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[3][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(3.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 770      //-----------------------------------------------------------------------
 771      //real part
 772      fDiffFlowGenFunRe4->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[4][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(4.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 773      //imaginary part
 774      fDiffFlowGenFunIm4->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[4][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(4.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 775      //-----------------------------------------------------------------------
 776      //real part
 777      fDiffFlowGenFunRe5->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[5][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(5.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 778      //imaginary part
 779      fDiffFlowGenFunIm5->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[5][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(5.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 780      //-----------------------------------------------------------------------
 781      //real part
 782      fDiffFlowGenFunRe6->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[6][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(6.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 783      //imaginary part
 784      fDiffFlowGenFunIm6->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[6][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(6.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 785      //-----------------------------------------------------------------------
 786      //real part
 787      fDiffFlowGenFunRe7->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[7][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(7.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 788      //imaginary part
 789      fDiffFlowGenFunIm7->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[7][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(7.+1.)/nSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.);
 790    }
 791    //}
 792   }
 793  }//ending the second loop over data
 794  */
 795
 796  }//end of if(nEventNSelTracksIntFlow>9)
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810  //off the record: numerical equations for cumulants solved up to different highest order
 811  if(fOtherEquations)
 812  {
 813   //running over data
 814   Int_t nPrimOE = anEvent->NumberOfTracks();//total multiplicity
 815
 816   Int_t nEventNSelTracksIntFlowOE = anEvent->GetEventNSelTracksIntFlow();
 817
 818   Double_t genfunG4[fgkPmax4][fgkQmax4];
 819   Double_t genfunG6[fgkPmax6][fgkQmax6];
 820   Double_t genfunG8[fgkPmax8][fgkQmax8];
 821   Double_t genfunG16[fgkPmax16][fgkQmax16];
 822   for(Int_t p=0;p<fgkPmax16;p++)
 823   {
 824    for(Int_t q=0;q<fgkQmax16;q++)
 825    {
 826     genfunG16[p][q]=1.;
 827     if(p<fgkPmax8 && q<fgkQmax8)
 828     {
 829      genfunG8[p][q]=1.;
 830      if(p<fgkPmax6 && q<fgkQmax6)
 831      {
 832       genfunG6[p][q]=1.;
 833       if(p<fgkPmax4 && q<fgkQmax4)
 834       {
 835        genfunG4[p][q]=1.;
 836       }
 837      }
 838     }
 839    }
 840   }
 841
 842   //multiplicities:
 843   if(nEventNSelTracksIntFlowOE>15) fAvMultIntFlow16GFC->Fill(0.,nEventNSelTracksIntFlowOE,1);
 844   if(nEventNSelTracksIntFlowOE>7) fAvMultIntFlow8GFC->Fill(0.,nEventNSelTracksIntFlowOE,1);
 845   if(nEventNSelTracksIntFlowOE>5) fAvMultIntFlow6GFC->Fill(0.,nEventNSelTracksIntFlowOE,1);
 846   if(nEventNSelTracksIntFlowOE>3) fAvMultIntFlow4GFC->Fill(0.,nEventNSelTracksIntFlowOE,1);
 847
 848   //first loop over data: evaluating the generating function G[p][q] for integrated flow
 849   for(Int_t i=0;i<nPrimOE;i++)
 850   {
 851    fTrack=anEvent->GetTrack(i);
 852    if(fTrack && fTrack->UseForIntegratedFlow())
 853    {
 854     for(Int_t p=0;p<fgkPmax16;p++)
 855     {
 856      for(Int_t q=0;q<fgkQmax16;q++)
 857      {
 858       if(nEventNSelTracksIntFlowOE>15)
 859       {
 860        genfunG16[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksIntFlowOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax16));
 861       }
 862       if(p<fgkPmax8 && q<fgkQmax8)
 863       {
 864        if(nEventNSelTracksIntFlowOE>7)
 865        {
 866         genfunG8[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksIntFlowOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax8));
 867        }
 868        if(p<fgkPmax6 && q<fgkQmax6)
 869        {
 870         if(nEventNSelTracksIntFlowOE>5)
 871         {
 872          genfunG6[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksIntFlowOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax6));
 873         }
 874         if(p<fgkPmax4 && q<fgkQmax4)
 875         {
 876          if(nEventNSelTracksIntFlowOE>3)
 877          {
 878           genfunG4[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksIntFlowOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax4));
 879          }
 880         }
 881        }
 882       }
 883      }
 884     }
 885    }//end of if(fTrack && fTrack->UseForIntegratedFlow())
 886   }//ending the loop over data
 887
 888  //storing the value of G[p][q] in 2D profile in order to get automatically the avarage <G[p][q]>
 889  for(Int_t p=0;p<fgkPmax16;p++)
 890  {
 891   for(Int_t q=0;q<fgkQmax16;q++)
 892   {
 893    if(nEventNSelTracksIntFlowOE>15) fIntFlowGenFun16->Fill((Double_t)p,(Double_t)q,genfunG16[p][q],1);
 894    if(p<fgkPmax8 && q<fgkQmax8)
 895    {
 896     if(nEventNSelTracksIntFlowOE>7) fIntFlowGenFun8->Fill((Double_t)p,(Double_t)q,genfunG8[p][q],1);
 897     if(p<fgkPmax6 && q<fgkQmax6)
 898     {
 899      if(nEventNSelTracksIntFlowOE>5) fIntFlowGenFun6->Fill((Double_t)p,(Double_t)q,genfunG6[p][q],1);
 900      if(p<fgkPmax4 && q<fgkQmax4)
 901      {
 902       if(nEventNSelTracksIntFlowOE>3) fIntFlowGenFun4->Fill((Double_t)p,(Double_t)q,genfunG4[p][q],1);
 903      }
 904     }
 905    }
 906   }
 907  }
 908 }//end of if(fOtherEquations)
 909
 910 }//end of Make()
 911
 912 //================================================================================================================
 913
 914 void AliFlowAnalysisWithCumulants::Finish()
 915 {
 916  //calculate the final results
 917  //AliCumulantsFunctions finalResults(fIntFlowGenFun,NULL,NULL, fIntFlowResults,fDiffFlowResults2,fDiffFlowResults4,fDiffFlowResults6,fDiffFlowResults8,fAvMultIntFlow,fQVectorComponents,  fQDist,fDiffFlowGenFunRe0,fDiffFlowGenFunRe1,fDiffFlowGenFunRe2, fDiffFlowGenFunRe3,fDiffFlowGenFunRe4,fDiffFlowGenFunRe5,fDiffFlowGenFunRe6,fDiffFlowGenFunRe7,fDiffFlowGenFunIm0,fDiffFlowGenFunIm1, fDiffFlowGenFunIm2,fDiffFlowGenFunIm3,fDiffFlowGenFunIm4,fDiffFlowGenFunIm5,fDiffFlowGenFunIm6,fDiffFlowGenFunIm7);
 918
 919  AliCumulantsFunctions finalResults(fIntFlowGenFun,fIntFlowGenFun4,fIntFlowGenFun6,fIntFlowGenFun8,fIntFlowGenFun16,fAvMultIntFlow4GFC, fAvMultIntFlow6GFC,fAvMultIntFlow8GFC,fAvMultIntFlow16GFC,fDiffFlowPtRPGenFunRe,fDiffFlowPtRPGenFunIm,fPtBinRPNoOfParticles, fDiffFlowEtaRPGenFunRe,fDiffFlowEtaRPGenFunIm,fEtaBinRPNoOfParticles,fDiffFlowPtPOIGenFunRe,fDiffFlowPtPOIGenFunIm,fPtBinPOINoOfParticles, fDiffFlowEtaPOIGenFunRe,fDiffFlowEtaPOIGenFunIm,fEtaBinPOINoOfParticles,fIntFlowResultsGFC,fDiffFlowResults2ndOrderGFC,fDiffFlowResults4thOrderGFC,fDiffFlowResults6thOrderGFC,fDiffFlowResults8thOrderGFC, fAvMultIntFlowGFC,fQVectorComponentsGFC,fAverageOfSquaredWeight,fCommonHistsResults2nd, fCommonHistsResults4th,fCommonHistsResults6th,fCommonHistsResults8th,fCommonHists);
 920
 921  finalResults.Calculate();
 922 }
 923
 924 //================================================================================================================
 925
 926 void AliFlowAnalysisWithCumulants::WriteHistograms(TString* outputFileName)
 927 {
 928  //store the final results in output .root file
 929  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 930  output->WriteObject(fHistList, "cobjGFC","SingleKey");
 931  delete output;
 932 }
 933
 934 //================================================================================================================
 935
 936 //================================================================================================================
 937
 938 void AliFlowAnalysisWithCumulants::WriteHistograms(TString outputFileName)
 939 {
 940  //store the final results in output .root file
 941  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 942  output->WriteObject(fHistList, "cobjGFC","SingleKey");
 943  delete output;
 944 }
 945
 946 //================================================================================================================
 947
 948
 949
 950