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