238f6dc421d18775d11337df03b66f01fe881f01
[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 {
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