uncertainties back in, incase of weights still behaving funny
[u/mrichter/AliRoot.git] / 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