f7afc6b2ab67e962bccebda63b27f036d1aaf31d
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithQCumulants.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 Q-cumulants * 
18  *                                * 
19  * author:  Ante Bilandzic        * 
20  *           (anteb@nikhef.nl)    *
21  *********************************/ 
22
23 #define AliFlowAnalysisWithQCumulants_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"
32 #include "TGraph.h"
33 #include "TParticle.h"
34 #include "TRandom3.h"
35 #include "TStyle.h"
36 #include "TProfile.h"
37 #include "TProfile2D.h" 
38 #include "TProfile3D.h"
39 #include "TMath.h"
40 #include "TArrow.h"
41 #include "TPaveLabel.h"
42 #include "TCanvas.h"
43 #include "AliFlowEventSimple.h"
44 #include "AliFlowTrackSimple.h"
45 #include "AliFlowAnalysisWithQCumulants.h"
46 #include "TArrayD.h"
47 #include "TRandom.h"
48
49 class TH1;
50 class TGraph;
51 class TPave;
52 class TLatex;
53 class TMarker;
54 class TRandom3;
55 class TObjArray;
56 class TList;
57 class TCanvas;
58 class TSystem;
59 class TROOT;
60 class AliFlowVector;
61 class TVector;
62
63 //================================================================================================================
64
65 ClassImp(AliFlowAnalysisWithQCumulants)
66
67 AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants():  
68  fTrack(NULL),
69  fHistList(NULL),
70  fDiffFlowList(NULL),
71  fWeightsList(NULL),
72  fAvMultIntFlowQC(NULL),
73  fQvectorComponents(NULL),
74  fIntFlowResultsQC(NULL),
75  fDiffFlowResults2ndOrderQC(NULL),
76  fDiffFlowResults4thOrderQC(NULL),
77  fCovariances(NULL),
78  fQvectorForEachEventX(NULL),//to be removed
79  fQvectorForEachEventY(NULL),//to be removed
80  fQCorrelations(NULL),
81  fWeightedQCorrelations(NULL),
82  fQProduct(NULL),
83  fDirectCorrelations(NULL),
84  f2PerPtBin1n1nPOI(NULL),
85  f4PerPtBin1n1n1n1nPOI(NULL),
86  f2PerEtaBin1n1nPOI(NULL),
87  f4PerEtaBin1n1n1n1nPOI(NULL), 
88  f2WPerPtBin1n1nPOI(NULL),
89  f4WPerPtBin1n1n1n1nPOI(NULL),
90  f2WPerEtaBin1n1nPOI(NULL),
91  f4WPerEtaBin1n1n1n1nPOI(NULL),
92  f2PerPtBin1n1nRP(NULL),
93  f4PerPtBin1n1n1n1nRP(NULL),
94  f2PerEtaBin1n1nRP(NULL),
95  f4PerEtaBin1n1n1n1nRP(NULL),
96  f2WPerPtBin1n1nRP(NULL),
97  f4WPerPtBin1n1n1n1nRP(NULL),
98  f2WPerEtaBin1n1nRP(NULL),
99  f4WPerEtaBin1n1n1n1nRP(NULL),
100  fCommonHists2nd(NULL),
101  fCommonHists4th(NULL),
102  fCommonHists6th(NULL),
103  fCommonHists8th(NULL),
104  fCommonHistsResults2nd(NULL),
105  fCommonHistsResults4th(NULL),
106  fCommonHistsResults6th(NULL),
107  fCommonHistsResults8th(NULL),
108  f2pDistribution(NULL),
109  f4pDistribution(NULL),
110  f6pDistribution(NULL),
111  f8pDistribution(NULL),
112  fnBinsPt(0),
113  fPtMin(0),
114  fPtMax(0),
115  fnBinsEta(0),
116  fEtaMin(0),
117  fEtaMax(0),
118  fEventCounter(0),
119  fUsePhiWeights(kFALSE),
120  fUsePtWeights(kFALSE),
121  fUseEtaWeights(kFALSE)
122 {
123  //constructor 
124  fHistList = new TList();
125  fDiffFlowList = new TList(); 
126  fDiffFlowList->SetName("DifferentialFlow"); 
127  fWeightsList = new TList();
128  fWeightsList->SetName("Weights");
129   
130  fnBinsPt = AliFlowCommonConstants::GetNbinsPt();
131  fPtMin   = AliFlowCommonConstants::GetPtMin();      
132  fPtMax   = AliFlowCommonConstants::GetPtMax();
133  
134  fnBinsEta = AliFlowCommonConstants::GetNbinsEta();
135  fEtaMin   = AliFlowCommonConstants::GetEtaMin();            
136  fEtaMax   = AliFlowCommonConstants::GetEtaMax();
137 }
138
139 AliFlowAnalysisWithQCumulants::~AliFlowAnalysisWithQCumulants()
140 {
141  //destructor
142  delete fHistList; 
143  delete fDiffFlowList;
144  delete fWeightsList; 
145 }
146
147 //================================================================================================================
148
149 void AliFlowAnalysisWithQCumulants::Init()
150 {
151  //various output histograms
152  //avarage multiplicity 
153  fAvMultIntFlowQC = new TProfile("fAvMultIntFlowQC","Average Multiplicity",1,0,1,"s");
154  fAvMultIntFlowQC->SetXTitle("");
155  fAvMultIntFlowQC->SetYTitle("");
156  fAvMultIntFlowQC->SetLabelSize(0.06);
157  fAvMultIntFlowQC->SetMarkerStyle(25);
158  fAvMultIntFlowQC->SetLabelOffset(0.01);
159  (fAvMultIntFlowQC->GetXaxis())->SetBinLabel(1,"Average Multiplicity");
160  fHistList->Add(fAvMultIntFlowQC);
161  
162  //Q-vector stuff
163  fQvectorComponents = new TProfile("fQvectorComponents","Avarage of Q-vector components",44,0.,44.,"s");
164  fQvectorComponents->SetXTitle("");
165  fQvectorComponents->SetYTitle("");
166  //fHistList->Add(fQvectorComponents);
167  
168  //final results for integrated flow from Q-cumulants
169  fIntFlowResultsQC = new TH1D("fIntFlowResultsQC","Integrated Flow from Q-cumulants",4,0,4);
170  //fIntFlowResults->SetXTitle("");
171  //fIntFlowResultsQC->SetYTitle("IntegFALSrated Flow");
172  fIntFlowResultsQC->SetLabelSize(0.06);
173  //fIntFlowResultsQC->SetTickLength(1);
174  fIntFlowResultsQC->SetMarkerStyle(25);
175  (fIntFlowResultsQC->GetXaxis())->SetBinLabel(1,"v_{n}{2}");
176  (fIntFlowResultsQC->GetXaxis())->SetBinLabel(2,"v_{n}{4}");
177  (fIntFlowResultsQC->GetXaxis())->SetBinLabel(3,"v_{n}{6}");
178  (fIntFlowResultsQC->GetXaxis())->SetBinLabel(4,"v_{n}{8}");
179  fHistList->Add(fIntFlowResultsQC);
180
181  //final results for differential flow from 2nd order Q-cumulant
182  fDiffFlowResults2ndOrderQC = new TH1D("fDiffFlowResults2ndOrderQC","Differential Flow from 2nd Order Q-cumulant",fnBinsPt,fPtMin,fPtMax);
183  fDiffFlowResults2ndOrderQC->SetXTitle("p_{t} [GeV]");
184  //fDiffFlowResults2ndOrderQC->SetYTitle("Differential Flow");
185  fHistList->Add(fDiffFlowResults2ndOrderQC);
186  
187  //final results for differential flow from 4th order Q-cumulant
188  fDiffFlowResults4thOrderQC = new TH1D("fDiffFlowResults4thOrderQC","Differential Flow from 4th Order Q-cumulant",fnBinsPt,fPtMin,fPtMax);
189  fDiffFlowResults4thOrderQC->SetXTitle("p_{t} [GeV]");
190  //fDiffFlowResults4thOrderQC->SetYTitle("Differential Flow");
191  fHistList->Add(fDiffFlowResults4thOrderQC);
192  
193  //final results for covariances (1st bin: <2*4>-<2>*<4>, 2nd bin: <2*6>-<2>*<6>, ...)
194  fCovariances = new TH1D("fCovariances","Covariances",6,0,6);
195  //fCovariances->SetXTitle("");
196  //fCovariances->SetYTitle("<covariance>");
197  fCovariances->SetLabelSize(0.04);
198  fCovariances->SetTickLength(1);
199  fCovariances->SetMarkerStyle(25);
200  (fCovariances->GetXaxis())->SetBinLabel(1,"Cov(2,4)");
201  (fCovariances->GetXaxis())->SetBinLabel(2,"Cov(2,6)");
202  (fCovariances->GetXaxis())->SetBinLabel(3,"Cov(2,8)");
203  (fCovariances->GetXaxis())->SetBinLabel(4,"Cov(4,6)");
204  (fCovariances->GetXaxis())->SetBinLabel(5,"Cov(4,8)");
205  (fCovariances->GetXaxis())->SetBinLabel(6,"Cov(6,8)");
206  fHistList->Add(fCovariances);
207   
208  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
209  //                   !!!! to be removed !!!! 
210  //profile containing the x-components of Q-vectors from all events 
211  fQvectorForEachEventX = new TProfile("fQvectorForEachEventX","x-components of Q-vectors",44000,1,44000,"s");
212  fHistList->Add(fQvectorForEachEventX);
213  
214  //profile containing the y-components of Q-vectors from all events 
215  fQvectorForEachEventY = new TProfile("fQvectorForEachEventY","y-components of Q-vectors",44000,1,44000,"s");
216  fHistList->Add(fQvectorForEachEventY);
217  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
218     
219  //multi-particle correlations calculated from Q-vectors
220  fQCorrelations = new TProfile("fQCorrelations","multi-particle correlations from Q-vectors",32,0,32,"s");
221  //fQCorrelations->SetXTitle("correlations");
222  //fQCorrelations->SetYTitle("");
223  fQCorrelations->SetTickLength(-0.01,"Y");
224  fQCorrelations->SetMarkerStyle(25);
225  fQCorrelations->SetLabelSize(0.03);
226  fQCorrelations->SetLabelOffset(0.01,"Y");
227  
228  (fQCorrelations->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
229  (fQCorrelations->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
230  (fQCorrelations->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
231  (fQCorrelations->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
232  
233  (fQCorrelations->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
234  (fQCorrelations->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
235  (fQCorrelations->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
236  (fQCorrelations->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
237  
238  (fQCorrelations->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
239  (fQCorrelations->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
240  (fQCorrelations->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
241  (fQCorrelations->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
242  (fQCorrelations->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
243  (fQCorrelations->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
244  (fQCorrelations->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
245  
246  (fQCorrelations->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
247  (fQCorrelations->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
248  (fQCorrelations->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
249  (fQCorrelations->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
250  
251  (fQCorrelations->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
252  (fQCorrelations->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
253  (fQCorrelations->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
254  (fQCorrelations->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
255  
256  (fQCorrelations->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
257
258  (fQCorrelations->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
259  
260  fHistList->Add(fQCorrelations);
261  
262  
263  
264  
265  //weighted multi-particle correlations calculated from Q-vectors
266  fWeightedQCorrelations = new TProfile("fWeightedQCorrelations","weighted multi-particle correlations from Q-vectors",100,0,100,"s");
267  //fWeightedQCorrelations->SetXTitle("correlations");
268  //fWeightedQCorrelations->SetYTitle("");
269  fWeightedQCorrelations->SetTickLength(-0.01,"Y");
270  fWeightedQCorrelations->SetMarkerStyle(25);
271  fWeightedQCorrelations->SetLabelSize(0.03);
272  fWeightedQCorrelations->SetLabelOffset(0.01,"Y");
273  
274  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(1,"<w_{1}w_{2}cos(n(#phi_{1}-#phi_{2}))>");
275  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(2,"<w_{1}^{2}w_{2}^{2}cos(2n(#phi_{1}-#phi_{2}))>");
276  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(3,"<w_{1}^{3}w_{2}^{3}cos(3n(#phi_{1}-#phi_{2}))>");
277  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(4,"<w_{1}^{4}w_{2}^{4}cos(4n(#phi_{1}-#phi_{2}))>");
278  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(5,"<w_{1}^{3}w_{2}cos(n(#phi_{1}-#phi_{2}))>");
279  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(6,"<w_{1}^{2}w_{2}w_{3}cos(n(#phi_{1}-#phi_{2}))>");
280  
281  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(11,"<w_{1}w_{2}w_{3}^{2}cos(n(2#phi_{1}-#phi_{2}-#phi_{3}))>");
282  
283  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(21,"<w_{1}w_{2}w_{3}w_{4}cos(n(#phi_{1}+#phi_{2}-#phi_{3}-#phi_{4}))>");
284  
285  /*
286  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
287  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
288  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
289  
290  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
291  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
292  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
293  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
294  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
295  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
296  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
297  
298  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
299  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
300  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
301  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
302  
303  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
304  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
305  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
306  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
307  
308  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
309
310  (fWeightedQCorrelations->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
311  */
312  
313  fHistList->Add(fWeightedQCorrelations);
314  
315  
316  
317  
318  //average products
319  fQProduct = new TProfile("fQProduct","average of products",6,0,6,"s");
320  fQProduct->SetTickLength(-0.01,"Y");
321  fQProduct->SetMarkerStyle(25);
322  fQProduct->SetLabelSize(0.03);
323  fQProduct->SetLabelOffset(0.01,"Y");
324  (fQProduct->GetXaxis())->SetBinLabel(1,"<<2*4>>");
325  (fQProduct->GetXaxis())->SetBinLabel(2,"<<2*6>>");
326  (fQProduct->GetXaxis())->SetBinLabel(3,"<<2*8>>");
327  (fQProduct->GetXaxis())->SetBinLabel(4,"<<4*6>>");
328  (fQProduct->GetXaxis())->SetBinLabel(5,"<<4*8>>");
329  (fQProduct->GetXaxis())->SetBinLabel(6,"<<6*8>>");
330  fQProduct->SetXTitle("");
331  fQProduct->SetYTitle("");
332  fHistList->Add(fQProduct);
333  
334  //weighted multi-particle correlations calculated with nested loops (0..100 integrated flow; 100..200 differential flow)
335  fDirectCorrelations = new TProfile("fDirectCorrelations","multi-particle correlations with nested loops",200,0,200,"s");
336  fDirectCorrelations->SetXTitle("");
337  fDirectCorrelations->SetYTitle("correlations");
338  fHistList->Add(fDirectCorrelations);
339  
340  //f2PerPtBin1n1nRP
341  f2PerPtBin1n1nRP = new TProfile("f2PerPtBin1n1nRP","<2'>_{n|n}",fnBinsPt,fPtMin,fPtMax,"s");
342  f2PerPtBin1n1nRP->SetXTitle("p_{t} [GeV]");
343  fDiffFlowList->Add(f2PerPtBin1n1nRP);
344  
345  //f4PerPtBin1n1n1n1nRP
346  f4PerPtBin1n1n1n1nRP = new TProfile("f4PerPtBin1n1n1n1nRP","<4'>_{n,n|n,n}",fnBinsPt,fPtMin,fPtMax,"s");
347  f4PerPtBin1n1n1n1nRP->SetXTitle("p_{t} [GeV]");
348  fDiffFlowList->Add(f4PerPtBin1n1n1n1nRP);
349  
350  //f2PerEtaBin1n1nRP
351  f2PerEtaBin1n1nRP = new TProfile("f2PerEtaBin1n1nRP","<2'>_{n|n}",fnBinsEta,fEtaMin,fEtaMax,"s");
352  f2PerEtaBin1n1nRP->SetXTitle("#eta");
353  fDiffFlowList->Add(f2PerEtaBin1n1nRP);
354  
355  //f4PerEtaBin1n1n1n1nRP
356  f4PerEtaBin1n1n1n1nRP = new TProfile("f4PerEtaBin1n1n1n1nRP","<4'>_{n,n|n,n}",fnBinsEta,fEtaMin,fEtaMax,"s");
357  f4PerEtaBin1n1n1n1nRP->SetXTitle("#eta");
358  fDiffFlowList->Add(f4PerEtaBin1n1n1n1nRP);
359  
360  //f2PerPtBin1n1nPOI
361  f2PerPtBin1n1nPOI = new TProfile("f2PerPtBin1n1nPOI","<2'>_{n|n}",fnBinsPt,fPtMin,fPtMax,"s");
362  f2PerPtBin1n1nPOI->SetXTitle("#eta");
363  fDiffFlowList->Add(f2PerPtBin1n1nPOI);
364  
365  //f4PerPtBin1n1n1n1nPOI
366  f4PerPtBin1n1n1n1nPOI = new TProfile("f4PerPtBin1n1n1n1nPOI","<4'>_{n,n|n,n}",fnBinsPt,fPtMin,fPtMax,"s");
367  f4PerPtBin1n1n1n1nPOI->SetXTitle("p_{t} [GeV]");
368  fDiffFlowList->Add(f4PerPtBin1n1n1n1nPOI);
369  
370  //f2PerEtaBin1n1nPOI
371  f2PerEtaBin1n1nPOI = new TProfile("f2PerEtaBin1n1nPOI","<2'>_{n|n}",fnBinsEta,fEtaMin,fEtaMax,"s");
372  f2PerEtaBin1n1nPOI->SetXTitle("#eta");
373  fDiffFlowList->Add(f2PerEtaBin1n1nPOI);
374  
375  //f4PerEtaBin1n1n1n1nPOI
376  f4PerEtaBin1n1n1n1nPOI = new TProfile("f4PerEtaBin1n1n1n1nPOI","<4'>_{n,n|n,n}",fnBinsEta,fEtaMin,fEtaMax,"s");
377  f4PerEtaBin1n1n1n1nPOI->SetXTitle("#eta");
378  fDiffFlowList->Add(f4PerEtaBin1n1n1n1nPOI);
379  
380  //f2WPerPtBin1n1nPOI
381  f2WPerPtBin1n1nPOI = new TProfile("f2WPerPtBin1n1nPOI","<2'>_{n|n}",fnBinsPt,fPtMin,fPtMax,"s");
382  f2WPerPtBin1n1nPOI->SetXTitle("#pt");
383  fDiffFlowList->Add(f2WPerPtBin1n1nPOI);
384  
385  //f4WPerPtBin1n1n1n1nPOI
386  f4WPerPtBin1n1n1n1nPOI = new TProfile("f4WPerPtBin1n1n1n1nPOI","<4'>_{n,n|n,n}",fnBinsPt,fPtMin,fPtMax,"s");
387  f4WPerPtBin1n1n1n1nPOI->SetXTitle("#Pt");
388  fDiffFlowList->Add(f4WPerPtBin1n1n1n1nPOI);
389  
390  //f2WPerEtaBin1n1nPOI
391  f2WPerEtaBin1n1nPOI = new TProfile("f2WPerEtaBin1n1nPOI","<2'>_{n|n}",fnBinsEta,fEtaMin,fEtaMax,"s");
392  f2WPerEtaBin1n1nPOI->SetXTitle("#eta");
393  fDiffFlowList->Add(f2WPerEtaBin1n1nPOI);
394  
395  //f4WPerEtaBin1n1n1n1nPOI
396  f4WPerEtaBin1n1n1n1nPOI = new TProfile("f4WPerEtaBin1n1n1n1nPOI","<4'>_{n,n|n,n}",fnBinsEta,fEtaMin,fEtaMax,"s");
397  f4WPerEtaBin1n1n1n1nPOI->SetXTitle("#eta");
398  fDiffFlowList->Add(f4WPerEtaBin1n1n1n1nPOI);
399  
400  //f2WPerPtBin1n1nRP
401  f2WPerPtBin1n1nRP = new TProfile("f2WPerPtBin1n1nRP","<2'>_{n|n}",fnBinsPt,fPtMin,fPtMax,"s");
402  f2WPerPtBin1n1nRP->SetXTitle("#pt");
403  fDiffFlowList->Add(f2WPerPtBin1n1nRP);
404  
405  //f4WPerPtBin1n1n1n1nRP
406  f4WPerPtBin1n1n1n1nRP = new TProfile("f4WPerPtBin1n1n1n1nRP","<4'>_{n,n|n,n}",fnBinsPt,fPtMin,fPtMax,"s");
407  f4WPerPtBin1n1n1n1nRP->SetXTitle("#Pt");
408  fDiffFlowList->Add(f4WPerPtBin1n1n1n1nRP);
409  
410  //f2WPerEtaBin1n1nRP
411  f2WPerEtaBin1n1nRP = new TProfile("f2WPerEtaBin1n1nRP","<2'>_{n|n}",fnBinsEta,fEtaMin,fEtaMax,"s");
412  f2WPerEtaBin1n1nRP->SetXTitle("#eta");
413  fDiffFlowList->Add(f2WPerEtaBin1n1nRP);
414  
415  //f4WPerEtaBin1n1n1n1nRP
416  f4WPerEtaBin1n1n1n1nRP = new TProfile("f4WPerEtaBin1n1n1n1nRP","<4'>_{n,n|n,n}",fnBinsEta,fEtaMin,fEtaMax,"s");
417  f4WPerEtaBin1n1n1n1nRP->SetXTitle("#eta");
418  fDiffFlowList->Add(f4WPerEtaBin1n1n1n1nRP);
419  
420  //common control histogram (2nd order)
421  fCommonHists2nd = new AliFlowCommonHist("AliFlowCommonHist2ndOrderQC");
422  fHistList->Add(fCommonHists2nd);  
423  
424  //common control histogram (4th order)
425  fCommonHists4th = new AliFlowCommonHist("AliFlowCommonHist4thOrderQC");
426  fHistList->Add(fCommonHists4th);  
427  
428  //common control histogram (6th order)
429  fCommonHists6th = new AliFlowCommonHist("AliFlowCommonHist6thOrderQC");
430  fHistList->Add(fCommonHists6th);  
431  
432  //common control histogram (8th order)
433  fCommonHists8th = new AliFlowCommonHist("AliFlowCommonHist8thOrderQC");
434  fHistList->Add(fCommonHists8th);  
435   
436  //common histograms for final results (2nd order)
437  fCommonHistsResults2nd = new AliFlowCommonHistResults("AliFlowCommonHistResults2ndOrderQC");
438  fHistList->Add(fCommonHistsResults2nd); 
439  
440  //common histograms for final results (4th order)
441  fCommonHistsResults4th = new AliFlowCommonHistResults("AliFlowCommonHistResults4thOrderQC");
442  fHistList->Add(fCommonHistsResults4th);
443  
444  //common histograms for final results (6th order)
445  fCommonHistsResults6th = new AliFlowCommonHistResults("AliFlowCommonHistResults6thOrderQC");
446  fHistList->Add(fCommonHistsResults6th); 
447  
448  //common histograms for final results (8th order)
449  fCommonHistsResults8th = new AliFlowCommonHistResults("AliFlowCommonHistResults8thOrderQC");
450  fHistList->Add(fCommonHistsResults8th); 
451  
452  //weighted <2>_{n|n} distribution
453  f2pDistribution = new TH1D("f2pDistribution","<2>_{n|n} distribution",100000,-0.02,0.1);
454  f2pDistribution->SetXTitle("<2>_{n|n}");
455  f2pDistribution->SetYTitle("Counts");
456  fHistList->Add(f2pDistribution);
457
458  //weighted <4>_{n,n|n,n} distribution
459  f4pDistribution = new TH1D("f4pDistribution","<4>_{n,n|n,n} distribution",100000,-0.00025,0.002);
460  f4pDistribution->SetXTitle("<4>_{n,n|n,n}");
461  f4pDistribution->SetYTitle("Counts");
462  fHistList->Add(f4pDistribution); 
463  
464  //weighted <6>_{n,n,n|n,n,n} distribution
465  f6pDistribution = new TH1D("f6pDistribution","<6>_{n,n,n|n,n,n} distribution",100000,-0.000005,0.000025);
466  f6pDistribution->SetXTitle("<6>_{n,n,n|n,n,n}");
467  f6pDistribution->SetYTitle("Counts");
468  fHistList->Add(f6pDistribution);
469  
470  //weighted <8>_{n,n,n,n|n,n,n,n} distribution
471  f8pDistribution = new TH1D("f8pDistribution","<8>_{n,n,n,n|n,n,n,n} distribution",100000,-0.000000001,0.00000001);
472  f8pDistribution->SetXTitle("<8>_{n,n,n,n|n,n,n,n}");
473  f8pDistribution->SetYTitle("Counts");
474  fHistList->Add(f8pDistribution);
475  
476  // add list fWeightsList with weights to the main list
477  fHistList->Add(fWeightsList);
478   
479  // add list fDiffFlowList with histograms and profiles needed for differential flow to the main list 
480  fHistList->Add(fDiffFlowList); 
481 }//end of Init()
482
483 //================================================================================================================
484
485 void AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
486 {
487  // running over data 
488  
489  Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = nRP + nPOI + rest  
490  Int_t nRP = anEvent->GetEventNSelTracksIntFlow(); // nRP = number of particles used to determine the reaction plane
491  
492  Int_t n = 2; // int flow harmonic (to be improved)
493  
494  //needed for debugging: (to be improved - add explanation here) 
495  //Bool_t bNestedLoops=kTRUE;
496  //if(!(bNestedLoops)||(nPrim>0&&nPrim<12))
497  //{
498  //if(nPrim>0&&nPrim<10)
499  //{
500  
501  
502  
503  //---------------------------------------------------------------------------------------------------------
504  // non-weighted and weighted Q-vectors of an event built-up from RP particles evaluated in harmonics n, 2n, 3n and 4n:
505  AliFlowVector afvQvector1n, afvQvector2n, afvQvector3n, afvQvector4n;
506
507  // non-weighted Q-vector in harmonic n: 
508  afvQvector1n.Set(0.,0.);
509  afvQvector1n.SetMult(0);
510  afvQvector1n = anEvent->GetQ(1*n); 
511  
512  // non-weighted Q-vector in harmonic 2n: 
513  afvQvector2n.Set(0.,0.);
514  afvQvector2n.SetMult(0);
515  afvQvector2n = anEvent->GetQ(2*n); // to be improved: weights   
516           
517  // non-weighted Q-vector in harmonic 3n:                                                                 
518  afvQvector3n.Set(0.,0.);
519  afvQvector3n.SetMult(0);
520  afvQvector3n = anEvent->GetQ(3*n); // to be improved: weights
521  
522  // non-weighted Q-vector in harmonic 4n:
523  afvQvector4n.Set(0.,0.);
524  afvQvector4n.SetMult(0);
525  afvQvector4n = anEvent->GetQ(4*n); // to be improved: weights
526             
527             
528  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
529  //                        !!!! to be removed !!!!
530  fQvectorForEachEventX->Fill(1.*(++fEventCounter),afvQvector1n.X());
531  fQvectorForEachEventY->Fill(1.*(fEventCounter),afvQvector1n.Y()); 
532  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
533                        
534                                   
535                                                         
536  //---------------------------------------------------------------------------------------------------------
537  
538  //multiplicity of RP particles:
539  Double_t dMult = afvQvector1n.GetMult(); // to be improved (name, this is actually weighted multiplicity)
540  
541  fAvMultIntFlowQC->Fill(0.,dMult,1.); // to be removed (this info is also stored in one of control histograms)
542  
543  //---------------------------------------------------------------------------------------------------------
544  //
545  //                                          *******************
546  //                                          **** Q-vectors ****
547  //                                          *******************
548  //
549  Double_t reQ2nQ1nstarQ1nstar = pow(afvQvector1n.X(),2.)*afvQvector2n.X()+2.*afvQvector1n.X()*afvQvector1n.Y()*afvQvector2n.Y()-pow(afvQvector1n.Y(),2.)*afvQvector2n.X();//Re[Q_{2n} Q_{n}^* Q_{n}^*]
550  //Double_t imQ2nQ1nstarQ1nstar = pow(Qvector1n.X(),2.)*Qvector2n.Y()-2.*Qvector1n.X()*Qvector1n.Y()*Qvector2n.X()-pow(Qvector1n.Y(),2.)*Qvector2n.Y();//Im[Q_{2n} Q_{n}^* Q_{n}^*]
551  Double_t reQ1nQ1nQ2nstar = reQ2nQ1nstarQ1nstar;//Re[Q_{n} Q_{n} Q_{2n}^*] = Re[Q_{2n} Q_{n}^* Q_{n}^*]
552  Double_t reQ3nQ1nQ2nstarQ2nstar = (pow(afvQvector2n.X(),2.)-pow(afvQvector2n.Y(),2.))*(afvQvector3n.X()*afvQvector1n.X()-afvQvector3n.Y()*afvQvector1n.Y())+2.*afvQvector2n.X()*afvQvector2n.Y()*(afvQvector3n.X()*afvQvector1n.Y()+afvQvector3n.Y()*afvQvector1n.X());
553  //Double_t imQ3nQ1nQ2nstarQ2nstar = calculate and implement this (deleteMe) 
554  Double_t reQ2nQ2nQ3nstarQ1nstar = reQ3nQ1nQ2nstarQ2nstar;
555  Double_t reQ4nQ2nstarQ2nstar = pow(afvQvector2n.X(),2.)*afvQvector4n.X()+2.*afvQvector2n.X()*afvQvector2n.Y()*afvQvector4n.Y()-pow(afvQvector2n.Y(),2.)*afvQvector4n.X();//Re[Q_{4n} Q_{2n}^* Q_{2n}^*]
556  //Double_t imQ4nQ2nstarQ2nstar = calculate and implement this (deleteMe)
557  Double_t reQ2nQ2nQ4nstar = reQ4nQ2nstarQ2nstar;
558  Double_t reQ4nQ3nstarQ1nstar = afvQvector4n.X()*(afvQvector3n.X()*afvQvector1n.X()-afvQvector3n.Y()*afvQvector1n.Y())+afvQvector4n.Y()*(afvQvector3n.X()*afvQvector1n.Y()+afvQvector3n.Y()*afvQvector1n.X());//Re[Q_{4n} Q_{3n}^* Q_{n}^*]
559  Double_t reQ3nQ1nQ4nstar = reQ4nQ3nstarQ1nstar;//Re[Q_{3n} Q_{n} Q_{4n}^*] = Re[Q_{4n} Q_{3n}^* Q_{n}^*]
560  //Double_t imQ4nQ3nstarQ1nstar = calculate and implement this (deleteMe)
561  Double_t reQ3nQ2nstarQ1nstar = afvQvector3n.X()*afvQvector2n.X()*afvQvector1n.X()-afvQvector3n.X()*afvQvector2n.Y()*afvQvector1n.Y()+afvQvector3n.Y()*afvQvector2n.X()*afvQvector1n.Y()+afvQvector3n.Y()*afvQvector2n.Y()*afvQvector1n.X();//Re[Q_{3n} Q_{2n}^* Q_{n}^*]
562  Double_t reQ2nQ1nQ3nstar = reQ3nQ2nstarQ1nstar;//Re[Q_{2n} Q_{n} Q_{3n}^*] = Re[Q_{3n} Q_{2n}^* Q_{n}^*]
563  //Double_t imQ3nQ2nstarQ1nstar; //calculate and implement this (deleteMe)
564  Double_t reQ3nQ1nstarQ1nstarQ1nstar = afvQvector3n.X()*pow(afvQvector1n.X(),3)-3.*afvQvector1n.X()*afvQvector3n.X()*pow(afvQvector1n.Y(),2)+3.*afvQvector1n.Y()*afvQvector3n.Y()*pow(afvQvector1n.X(),2)-afvQvector3n.Y()*pow(afvQvector1n.Y(),3);//Re[Q_{3n} Q_{n}^* Q_{n}^* Q_{n}^*]
565  //Double_t imQ3nQ1nstarQ1nstarQ1nstar; //calculate and implement this (deleteMe)
566  Double_t xQ2nQ1nQ2nstarQ1nstar = pow(afvQvector2n.Mod()*afvQvector1n.Mod(),2);//|Q_{2n}|^2 |Q_{n}|^2
567  Double_t reQ4nQ2nstarQ1nstarQ1nstar = (afvQvector4n.X()*afvQvector2n.X()+afvQvector4n.Y()*afvQvector2n.Y())*(pow(afvQvector1n.X(),2)-pow(afvQvector1n.Y(),2))+2.*afvQvector1n.X()*afvQvector1n.Y()*(afvQvector4n.Y()*afvQvector2n.X()-afvQvector4n.X()*afvQvector2n.Y());//Re[Q_{4n} Q_{2n}^* Q_{n}^* Q_{n}^*] 
568  //Double_t imQ4nQ2nstarQ1nstarQ1nstar; //calculate and implement this (deleteMe)
569  Double_t reQ2nQ1nQ1nstarQ1nstarQ1nstar = (afvQvector2n.X()*afvQvector1n.X()-afvQvector2n.Y()*afvQvector1n.Y())*(pow(afvQvector1n.X(),3)-3.*afvQvector1n.X()*pow(afvQvector1n.Y(),2))+(afvQvector2n.X()*afvQvector1n.Y()+afvQvector1n.X()*afvQvector2n.Y())*(3.*afvQvector1n.Y()*pow(afvQvector1n.X(),2)-pow(afvQvector1n.Y(),3));//Re[Q_{2n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^*]
570  //Double_t imQ2nQ1nQ1nstarQ1nstarQ1nstar; //calculate and implement this (deleteMe)
571  Double_t reQ2nQ2nQ2nstarQ1nstarQ1nstar = pow(afvQvector2n.Mod(),2.)*(afvQvector2n.X()*(pow(afvQvector1n.X(),2.)-pow(afvQvector1n.Y(),2.))+2.*afvQvector2n.Y()*afvQvector1n.X()*afvQvector1n.Y());//Re[Q_{2n} Q_{2n} Q_{2n}^* Q_{n}^* Q_{n}^*]
572  //Double_t imQ2nQ2nQ2nstarQ1nstarQ1nstar = pow(Qvector2n.Mod(),2.)*(Qvector2n.Y()*(pow(Qvector1n.X(),2.)-pow(Qvector1n.Y(),2.))-2.*Qvector2n.X()*Qvector1n.X()*Qvector1n.Y());//Im[Q_{2n} Q_{2n} Q_{2n}^* Q_{n}^* Q_{n}^*]
573  Double_t reQ4nQ1nstarQ1nstarQ1nstarQ1nstar = pow(afvQvector1n.X(),4.)*afvQvector4n.X()-6.*pow(afvQvector1n.X(),2.)*afvQvector4n.X()*pow(afvQvector1n.Y(),2.)+pow(afvQvector1n.Y(),4.)*afvQvector4n.X()+4.*pow(afvQvector1n.X(),3.)*afvQvector1n.Y()*afvQvector4n.Y()-4.*pow(afvQvector1n.Y(),3.)*afvQvector1n.X()*afvQvector4n.Y();//Re[Q_{4n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
574  //Double_t imQ4nQ1nstarQ1nstarQ1nstarQ1nstar = pow(Qvector1n.X(),4.)*Qvector4n.Y()-6.*pow(Qvector1n.X(),2.)*Qvector4n.Y()*pow(Qvector1n.Y(),2.)+pow(Qvector1n.Y(),4.)*Qvector4n.Y()+4.*pow(Qvector1n.Y(),3.)*Qvector1n.X()*Qvector4n.X()-4.*pow(Qvector1n.X(),3.)*Qvector1n.Y()*Qvector4n.X();//Im[Q_{4n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
575  Double_t reQ3nQ1nQ2nstarQ1nstarQ1nstar = pow(afvQvector1n.Mod(),2.)*(afvQvector1n.X()*afvQvector2n.X()*afvQvector3n.X()-afvQvector3n.X()*afvQvector1n.Y()*afvQvector2n.Y()+afvQvector2n.X()*afvQvector1n.Y()*afvQvector3n.Y()+afvQvector1n.X()*afvQvector2n.Y()*afvQvector3n.Y());//Re[Q_{3n} Q_{n} Q_{2n}^* Q_{n}^* Q_{n}^*]
576  //Double_t imQ3nQ1nQ2nstarQ1nstarQ1nstar = pow(afvQvector1n.Mod(),2.)*(-afvQvector2n.X()*afvQvector3n.X()*afvQvector1n.Y()-afvQvector1n.X()*afvQvector3n.X()*afvQvector2n.Y()+afvQvector1n.X()*afvQvector2n.X()*afvQvector3n.Y()-afvQvector1n.Y()*afvQvector2n.Y()*afvQvector3n.Y());//Im[Q_{3n} Q_{n} Q_{2n}^* Q_{n}^* Q_{n}^*]
577  Double_t reQ2nQ2nQ1nstarQ1nstarQ1nstarQ1nstar = (pow(afvQvector1n.X(),2.)*afvQvector2n.X()-2.*afvQvector1n.X()*afvQvector2n.X()*afvQvector1n.Y()-afvQvector2n.X()*pow(afvQvector1n.Y(),2.)+afvQvector2n.Y()*pow(afvQvector1n.X(),2.)+2.*afvQvector1n.X()*afvQvector1n.Y()*afvQvector2n.Y()-pow(afvQvector1n.Y(),2.)*afvQvector2n.Y())*(pow(afvQvector1n.X(),2.)*afvQvector2n.X()+2.*afvQvector1n.X()*afvQvector2n.X()*afvQvector1n.Y()-afvQvector2n.X()*pow(afvQvector1n.Y(),2.)-afvQvector2n.Y()*pow(afvQvector1n.X(),2.)+2.*afvQvector1n.X()*afvQvector1n.Y()*afvQvector2n.Y()+pow(afvQvector1n.Y(),2.)*afvQvector2n.Y());//Re[Q_{2n} Q_{2n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
578  //Double_t imQ2nQ2nQ1nstarQ1nstarQ1nstarQ1nstar = 2.*(pow(afvQvector1n.X(),2.)*afvQvector2n.X()-afvQvector2n.X()*pow(afvQvector1n.Y(),2.)+2.*afvQvector1n.X()*afvQvector1n.Y()*afvQvector2n.Y())*(pow(afvQvector1n.X(),2.)*afvQvector2n.Y()-2.*afvQvector1n.X()*afvQvector1n.Y()*afvQvector2n.X()-pow(afvQvector1n.Y(),2.)*afvQvector2n.Y());//Im[Q_{2n} Q_{2n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
579  Double_t reQ3nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = pow(afvQvector1n.Mod(),2.)*(pow(afvQvector1n.X(),3.)*afvQvector3n.X()-3.*afvQvector1n.X()*afvQvector3n.X()*pow(afvQvector1n.Y(),2.)+3.*pow(afvQvector1n.X(),2.)*afvQvector1n.Y()*afvQvector3n.Y()-pow(afvQvector1n.Y(),3.)*afvQvector3n.Y());//Re[Q_{3n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
580  //Double_t imQ3nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = pow(afvQvector1n.Mod(),2.)*(pow(afvQvector1n.Y(),3.)*afvQvector3n.X()-3.*afvQvector1n.Y()*afvQvector3n.X()*pow(afvQvector1n.X(),2.)-3.*pow(afvQvector1n.Y(),2.)*afvQvector1n.X()*afvQvector3n.Y()+pow(afvQvector1n.X(),3.)*afvQvector3n.Y());//Im[Q_{3n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
581  Double_t xQ2nQ1nQ1nQ2nstarQ1nstarQ1nstar = pow(afvQvector2n.Mod(),2.)*pow(afvQvector1n.Mod(),4.);//|Q_{2n}|^2 |Q_{n}|^4
582  Double_t reQ2nQ1nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = pow(afvQvector1n.Mod(),4.)*(pow(afvQvector1n.X(),2.)*afvQvector2n.X()-afvQvector2n.X()*pow(afvQvector1n.Y(),2.)+2.*afvQvector1n.X()*afvQvector1n.Y()*afvQvector2n.Y());//Re[Q_{2n} Q_{n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
583  //Double_t imQ2nQ1nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = pow(afvQvector1n.Mod(),4.)*(pow(afvQvector1n.X(),2.)*afvQvector2n.Y()-afvQvector2n.Y()*pow(afvQvector1n.Y(),2.)-2.*afvQvector1n.X()*afvQvector2n.X()*afvQvector1n.Y());//Re[Q_{2n} Q_{n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
584  //---------------------------------------------------------------------------------------------------------
585  
586  //---------------------------------------------------------------------------------------------------------
587  //
588  //                                        **************************************
589  //                                        **** multi-particle correlations: ****
590  //                                        **************************************
591  //
592  // Remark 1: multi-particle correlations calculated with Q-vectors are stored in fQCorrelations.
593  // Remark 2: binning of fQCorrelations is organized as follows: 
594  //
595  // 1st bin: <2>_{n|n} = two1n1n
596  // 2nd bin: <2>_{2n|2n} = two2n2n
597  // 3rd bin: <2>_{3n|3n} = two3n3n
598  // 4th bin: <2>_{4n|4n} = two4n4n
599  // 5th bin: --  EMPTY --
600  // 6th bin: <3>_{2n|n,n} = three2n1n1n
601  // 7th bin: <3>_{3n|2n,n} = three3n2n1n
602  // 8th bin: <3>_{4n|2n,2n} = three4n2n2n
603  // 9th bin: <3>_{4n|3n,n} = three4n3n1n
604  //10th bin: --  EMPTY --
605  //11th bin: <4>_{n,n|n,n} = four1n1n1n1n
606  //12th bin: <4>_{2n,n|2n,n} = four2n1n2n1n
607  //13th bin: <4>_{2n,2n|2n,2n} = four2n2n2n2n
608  //14th bin: <4>_{3n|n,n,n} = four3n1n1n1n
609  //15th bin: <4>_{3n,n|3n,n} = four3n1n3n1n 
610  //16th bin: <4>_{3n,n|2n,2n} = four3n1n2n2n
611  //17th bin: <4>_{4n|2n,n,n} = four4n2n1n1n
612  //18th bin: --  EMPTY --
613  //19th bin: <5>_{2n|n,n,n,n} = five2n1n1n1n1n
614  //20th bin: <5>_{2n,2n|2n,n,n} = five2n2n2n1n1n
615  //21st bin: <5>_{3n,n|2n,n,n} = five3n1n2n1n1n
616  //22nd bin: <5>_{4n|n,n,n,n} = five4n1n1n1n1n  
617  //23rd bin: --  EMPTY --
618  //24th bin: <6>_{n,n,n|n,n,n} = six1n1n1n1n1n1n
619  //25th bin: <6>_{2n,n,n|2n,n,n} = six2n1n1n2n1n1n
620  //26th bin: <6>_{2n,2n|n,n,n,n} = six2n2n1n1n1n1n
621  //27th bin: <6>_{3n,n|n,n,n,n} = six3n1n1n1n1n1n
622  //28th bin: --  EMPTY --
623  //29th bin: <7>_{2n,n,n|n,n,n,n} = seven2n1n1n1n1n1n1n
624  //30th bin: --  EMPTY --
625  //31st bin: <8>_{n,n,n,n|n,n,n,n} = eight1n1n1n1n1n1n1n1n
626
627  
628  // binning of fQProduct (all correlations are evaluated in harmonic n): 
629  // 1st bin: <2>*<4>
630  // 2nd bin: <2>*<6>
631  // 3rd bin: <2>*<8> 
632  // 4th bin: <4>*<6>
633  // 5th bin: <4>*<8>
634  // 6th bin: <6>*<8>
635          
636  // 2-particle
637  Double_t two1n1n = 0., two2n2n = 0., two3n3n = 0., two4n4n = 0.; 
638  
639  if(dMult>1)
640  {
641   //fill the common control histogram (2nd order): 
642   fCommonHists2nd->FillControlHistograms(anEvent); 
643   
644   two1n1n = (pow(afvQvector1n.Mod(),2.)-dMult)/(dMult*(dMult-1.)); // <2>_{n|n} = <cos(n*(phi1-phi2))>
645   two2n2n = (pow(afvQvector2n.Mod(),2.)-dMult)/(dMult*(dMult-1.)); // <2>_{2n|2n} = <cos(2n*(phi1-phi2))>
646   two3n3n = (pow(afvQvector3n.Mod(),2.)-dMult)/(dMult*(dMult-1.)); // <2>_{3n|3n} = <cos(3n*(phi1-phi2))>
647   two4n4n = (pow(afvQvector4n.Mod(),2.)-dMult)/(dMult*(dMult-1.)); // <2>_{4n|4n} = <cos(4n*(phi1-phi2))>
648     
649   fQCorrelations->Fill(0.,two1n1n,dMult*(dMult-1.));  
650   fQCorrelations->Fill(1.,two2n2n,dMult*(dMult-1.)); 
651   fQCorrelations->Fill(2.,two3n3n,dMult*(dMult-1.)); 
652   fQCorrelations->Fill(3.,two4n4n,dMult*(dMult-1.)); 
653   
654   f2pDistribution->Fill(two1n1n,dMult*(dMult-1.)); 
655  }
656  
657  // 3-particle
658  Double_t three2n1n1n=0., three3n2n1n=0., three4n2n2n=0., three4n3n1n=0.;
659  if(dMult>2)
660  {
661   three2n1n1n = (reQ2nQ1nstarQ1nstar-2.*pow(afvQvector1n.Mod(),2.)-pow(afvQvector2n.Mod(),2.)+2.*dMult)/(dMult*(dMult-1.)*(dMult-2.)); //Re[<3>_{2n|n,n}] = Re[<3>_{n,n|2n}] = <cos(n*(2.*phi1-phi2-phi3))>
662   three3n2n1n = (reQ3nQ2nstarQ1nstar-pow(afvQvector3n.Mod(),2.)-pow(afvQvector2n.Mod(),2.)-pow(afvQvector1n.Mod(),2.)+2.*dMult)/(dMult*(dMult-1.)*(dMult-2.)); //Re[<3>_{3n|2n,n}] = Re[<3>_{2n,n|3n}] = <cos(n*(3.*phi1-2.*phi2-phi3))>
663   three4n2n2n = (reQ4nQ2nstarQ2nstar-2.*pow(afvQvector2n.Mod(),2.)-pow(afvQvector4n.Mod(),2.)+2.*dMult)/(dMult*(dMult-1.)*(dMult-2.)); //Re[<3>_{4n|2n,2n}] = Re[<3>_{2n,2n|4n}] = <cos(n*(4.*phi1-2.*phi2-2.*phi3))>
664   three4n3n1n = (reQ4nQ3nstarQ1nstar-pow(afvQvector4n.Mod(),2.)-pow(afvQvector3n.Mod(),2.)-pow(afvQvector1n.Mod(),2.)+2.*dMult)/(dMult*(dMult-1.)*(dMult-2.)); //Re[<3>_{4n|3n,n}] = Re[<3>_{3n,n|4n}] = <cos(n*(4.*phi1-3.*phi2-phi3))>
665  
666   fQCorrelations->Fill(5.,three2n1n1n,dMult*(dMult-1.)*(dMult-2.)); 
667   fQCorrelations->Fill(6.,three3n2n1n,dMult*(dMult-1.)*(dMult-2.));
668   fQCorrelations->Fill(7.,three4n2n2n,dMult*(dMult-1.)*(dMult-2.)); 
669   fQCorrelations->Fill(8.,three4n3n1n,dMult*(dMult-1.)*(dMult-2.));    
670  }
671  
672  //4-particle
673  Double_t four1n1n1n1n=0., four2n2n2n2n=0., four2n1n2n1n=0., four3n1n1n1n=0., four4n2n1n1n=0., four3n1n2n2n=0., four3n1n3n1n=0.;  
674  if(dMult>3)
675  {
676   //fill the common control histogram (4th order):
677   fCommonHists4th->FillControlHistograms(anEvent); 
678
679   four1n1n1n1n = (2.*dMult*(dMult-3.)+pow(afvQvector1n.Mod(),4.)-4.*(dMult-2.)*pow(afvQvector1n.Mod(),2.)-2.*reQ2nQ1nstarQ1nstar+pow(afvQvector2n.Mod(),2.))/(dMult*(dMult-1)*(dMult-2.)*(dMult-3.));//<4>_{n,n|n,n}
680   four2n2n2n2n = (2.*dMult*(dMult-3.)+pow(afvQvector2n.Mod(),4.)-4.*(dMult-2.)*pow(afvQvector2n.Mod(),2.)-2.*reQ4nQ2nstarQ2nstar+pow(afvQvector4n.Mod(),2.))/(dMult*(dMult-1)*(dMult-2.)*(dMult-3.));//<4>_{2n,2n|2n,2n}
681   four2n1n2n1n = (xQ2nQ1nQ2nstarQ1nstar-2.*reQ3nQ2nstarQ1nstar-2.*reQ2nQ1nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))-((dMult-5.)*pow(afvQvector1n.Mod(),2.)+(dMult-4.)*pow(afvQvector2n.Mod(),2.)-pow(afvQvector3n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))+(dMult-6.)/((dMult-1.)*(dMult-2.)*(dMult-3.));//Re[<4>_{2n,n|2n,n}]
682   four3n1n1n1n = (reQ3nQ1nstarQ1nstarQ1nstar-3.*reQ3nQ2nstarQ1nstar-3.*reQ2nQ1nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))+(2.*pow(afvQvector3n.Mod(),2.)+3.*pow(afvQvector2n.Mod(),2.)+6.*pow(afvQvector1n.Mod(),2.)-6.*dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));//Re[<4>_{3n|n,n,n}]
683   four4n2n1n1n = (reQ4nQ2nstarQ1nstarQ1nstar-2.*reQ4nQ3nstarQ1nstar-reQ4nQ2nstarQ2nstar-2.*reQ3nQ2nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))-(reQ2nQ1nstarQ1nstar-2.*pow(afvQvector4n.Mod(),2.)-2.*pow(afvQvector3n.Mod(),2.)-3.*pow(afvQvector2n.Mod(),2.)-4.*pow(afvQvector1n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))-(6.)/((dMult-1.)*(dMult-2.)*(dMult-3.));//Re[<4>_{4n|2n,n,n}]
684   four3n1n2n2n = (reQ3nQ1nQ2nstarQ2nstar-reQ4nQ2nstarQ2nstar-reQ3nQ1nQ4nstar-2.*reQ3nQ2nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))-(2.*reQ1nQ1nQ2nstar-pow(afvQvector4n.Mod(),2.)-2.*pow(afvQvector3n.Mod(),2.)-4.*pow(afvQvector2n.Mod(),2.)-4.*pow(afvQvector1n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))-(6.)/((dMult-1.)*(dMult-2.)*(dMult-3.));//Re[<4>_{3n,n|2n,2n}] 
685   four3n1n3n1n = (pow(afvQvector3n.Mod(),2.)*pow(afvQvector1n.Mod(),2.)-2.*reQ4nQ3nstarQ1nstar-2.*reQ3nQ2nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))+(pow(afvQvector4n.Mod(),2.)-(dMult-4.)*pow(afvQvector3n.Mod(),2.)+pow(afvQvector2n.Mod(),2.)-(dMult-4.)*pow(afvQvector1n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))+(dMult-6.)/((dMult-1.)*(dMult-2.)*(dMult-3.));//<4>_{3n,n|3n,n}
686   //four_3n1n3n1n = Q3nQ1nQ3nstarQ1nstar/(M*(M-1.)*(M-2.)*(M-3.))-(2.*three_3n2n1n+2.*three_4n3n1n)/(M-3.)-(two_4n4n+M*two_3n3n+two_2n2n+M*two_1n1n)/((M-2.)*(M-3.))-M/((M-1.)*(M-2.)*(M-3.));//<4>_{3n,n|3n,n}
687   
688   fQCorrelations->Fill(10.,four1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
689   fQCorrelations->Fill(11.,four2n1n2n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
690   fQCorrelations->Fill(12.,four2n2n2n2n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
691   fQCorrelations->Fill(13.,four3n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
692   fQCorrelations->Fill(14.,four3n1n3n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
693   fQCorrelations->Fill(15.,four3n1n2n2n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));  
694   fQCorrelations->Fill(16.,four4n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)); 
695   
696   f4pDistribution->Fill(four1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
697   
698   fQProduct->Fill(0.,two1n1n*four1n1n1n1n,dMult*(dMult-1.)*dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
699  }
700
701  //5-particle
702  Double_t five2n1n1n1n1n=0., five2n2n2n1n1n=0., five3n1n2n1n1n=0., five4n1n1n1n1n=0.;
703  if(dMult>4)
704  {
705   five2n1n1n1n1n = (reQ2nQ1nQ1nstarQ1nstarQ1nstar-reQ3nQ1nstarQ1nstarQ1nstar+6.*reQ3nQ2nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-(reQ2nQ1nQ3nstar+3.*(dMult-6.)*reQ2nQ1nstarQ1nstar+3.*reQ1nQ1nQ2nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-(2.*pow(afvQvector3n.Mod(),2.)+3.*pow(afvQvector2n.Mod()*afvQvector1n.Mod(),2.)-3.*(dMult-4.)*pow(afvQvector2n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-3.*(pow(afvQvector1n.Mod(),4.)-2.*(2*dMult-5.)*pow(afvQvector1n.Mod(),2.)+2.*dMult*(dMult-4.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));//Re[<5>_{2n,n|n,n,n}]
706   
707   five2n2n2n1n1n = (reQ2nQ2nQ2nstarQ1nstarQ1nstar-reQ4nQ2nstarQ1nstarQ1nstar-2.*reQ2nQ2nQ3nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))+2.*(reQ4nQ2nstarQ2nstar+4.*reQ3nQ2nstarQ1nstar+reQ3nQ1nQ4nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))+(reQ2nQ2nQ4nstar-2.*(dMult-5.)*reQ2nQ1nstarQ1nstar+2.*reQ1nQ1nQ2nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-(2.*pow(afvQvector4n.Mod(),2.)+4.*pow(afvQvector3n.Mod(),2.)+1.*pow(afvQvector2n.Mod(),4.)-2.*(3.*dMult-10.)*pow(afvQvector2n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-(4.*pow(afvQvector1n.Mod(),2.)*pow(afvQvector2n.Mod(),2.)-4.*(dMult-5.)*pow(afvQvector1n.Mod(),2.)+4.*dMult*(dMult-6.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));//Re[<5>_{2n,2n|2n,n,n}]  
708
709   //five_2n2n2n1n1n = reQ2nQ2nQ2nstarQ1nstarQ1nstar/(M*(M-1.)*(M-2.)*(M-3.)*(M-4.))-(4.*four_2n1n2n1n+2.*four_3n1n2n2n+1.*four_2n2n2n2n+four_4n2n1n1n)/(M-4.)-(2.*three_4n3n1n+three_4n2n2n+three_4n2n2n+2.*three_3n2n1n)/((M-3.)*(M-4.))-(4.*three_3n2n1n+(2.*M-1.)*three_2n1n1n+2.*three_2n1n1n)/((M-3.)*(M-4.))-(two_4n4n+2.*two_3n3n+4.*(M-1.)*two_2n2n+2.*(2.*M-1.)*two_1n1n)/((M-2.)*(M-3.)*(M-4.))-(2.*M-1.)/((M-1.)*(M-2.)*(M-3.)*(M-4.)); //OK! 
710    
711   five4n1n1n1n1n = (reQ4nQ1nstarQ1nstarQ1nstarQ1nstar-6.*reQ4nQ2nstarQ1nstarQ1nstar-4.*reQ3nQ1nstarQ1nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))+(8.*reQ4nQ3nstarQ1nstar+3.*reQ4nQ2nstarQ2nstar+12.*reQ3nQ2nstarQ1nstar+12.*reQ2nQ1nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-(6.*pow(afvQvector4n.Mod(),2.)+8.*pow(afvQvector3n.Mod(),2.)+12.*pow(afvQvector2n.Mod(),2.)+24.*pow(afvQvector1n.Mod(),2.)-24.*dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));//Re[<5>_{4n|n,n,n,n}] 
712   
713   //five_4n1n1n1n1n = reQ4nQ1nstarQ1nstarQ1nstarQ1nstar/(M*(M-1.)*(M-2.)*(M-3.)*(M-4.)) -  (4.*four_3n1n1n1n+6.*four_4n2n1n1n)/(M-4.)  -  (6.*three_2n1n1n  + 12.*three_3n2n1n + 4.*three_4n3n1n + 3.*three_4n2n2n)/((M-3.)*(M-4.))  -  (4.*two_1n1n + 6.*two_2n2n + 4.*two_3n3n + 1.*two_4n4n)/((M-2.)*(M-3.)*(M-4.)) - 1./((M-1.)*(M-2.)*(M-3.)*(M-4.)); //OK!
714   
715   five3n1n2n1n1n = (reQ3nQ1nQ2nstarQ1nstarQ1nstar-reQ4nQ2nstarQ1nstarQ1nstar-reQ3nQ1nstarQ1nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-(reQ3nQ1nQ2nstarQ2nstar-3.*reQ4nQ3nstarQ1nstar-reQ4nQ2nstarQ2nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-((2.*dMult-13.)*reQ3nQ2nstarQ1nstar-reQ3nQ1nQ4nstar-9.*reQ2nQ1nstarQ1nstar)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-(2.*reQ1nQ1nQ2nstar+2.*pow(afvQvector4n.Mod(),2.)-2.*(dMult-5.)*pow(afvQvector3n.Mod(),2.)+2.*pow(afvQvector3n.Mod(),2.)*pow(afvQvector1n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))+(2.*(dMult-6.)*pow(afvQvector2n.Mod(),2.)-2.*pow(afvQvector2n.Mod(),2.)*pow(afvQvector1n.Mod(),2.)-pow(afvQvector1n.Mod(),4.)+2.*(3.*dMult-11.)*pow(afvQvector1n.Mod(),2.))/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))-4.*(dMult-6.)/((dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));//Re[<5>_{3n,n|2n,n,n}] 
716   
717   //five3n1n2n1n1n = reQ3nQ1nQ2nstarQ1nstarQ1nstar/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)) -  (four4n2n1n1n+four1n1n1n1n+four3n1n1n1n+2.*four2n1n2n1n+2.*three3n2n1n+2.*four3n1n3n1n+four3n1n2n2n)/(dMult-4.)  -   (2.*three4n3n1n+three4n2n2n+6.*three3n2n1n+three4n3n1n+2.*three3n2n1n+3.*three2n1n1n+2.*three2n1n1n+4.*two1n1n+2.*two2n2n+2.*two3n3n)/((dMult-3.)*(dMult-4.))  -  (5.*two1n1n + 4.*two2n2n + 3.*two3n3n + 1.*two4n4n + 2.)/((dMult-2.)*(dMult-3.)*(dMult-4.))  - 1./((dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)) ;//Re[<5>_{3n,n|2n,n,n}] //OK!
718   
719   //five3n1n2n1n1n = reQ3nQ1nQ2nstarQ1nstarQ1nstar/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)) -  (four4n2n1n1n+four1n1n1n1n+four3n1n1n1n+2.*four2n1n2n1n+2.*four3n1n3n1n+four3n1n2n2n)/(dMult-4.)  -      (2.*three4n3n1n+three4n2n2n+2.*dMult*three3n2n1n+three4n3n1n+2.*three3n2n1n+3.*three2n1n1n+2.*three2n1n1n)/((dMult-3.)*(dMult-4.))  -  ((4.*dMult-3.)*two1n1n + 2.*dMult*two2n2n + (2.*dMult-1.)*two3n3n + two4n4n)/((dMult-2.)*(dMult-3.)*(dMult-4.))  - (2.*dMult-1.)/((dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)) ;//Re[<5>_{3n,n|2n,n,n}] //OK!
720    
721   fQCorrelations->Fill(18.,five2n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)); 
722   fQCorrelations->Fill(19.,five2n2n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
723   fQCorrelations->Fill(20.,five3n1n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
724   fQCorrelations->Fill(21.,five4n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
725  }
726
727  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
728  //      !!!! to be removed: temporary fix !!!!
729  if(dMult>1)
730  {
731   two1n1n = (pow(afvQvector1n.Mod(),2.)-dMult)/(dMult*(dMult-1.));
732  }  
733  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
734     
735  //6-particle
736  Double_t six1n1n1n1n1n1n=0., six2n2n1n1n1n1n=0., six3n1n1n1n1n1n=0., six2n1n1n2n1n1n=0.;
737  if(dMult>5)
738  {
739   //fill the common control histogram (6th order):
740   fCommonHists6th->FillControlHistograms(anEvent); 
741  
742   six1n1n1n1n1n1n = (pow(afvQvector1n.Mod(),6.)+9.*xQ2nQ1nQ2nstarQ1nstar-6.*reQ2nQ1nQ1nstarQ1nstarQ1nstar)/(dMult*(dMult-1)*(dMult-2)*(dMult-3)*(dMult-4)*(dMult-5))+4.*(reQ3nQ1nstarQ1nstarQ1nstar-3.*reQ3nQ2nstarQ1nstar)/(dMult*(dMult-1)*(dMult-2)*(dMult-3)*(dMult-4)*(dMult-5))+2.*(9.*(dMult-4.)*reQ2nQ1nstarQ1nstar+2.*pow(afvQvector3n.Mod(),2.))/(dMult*(dMult-1)*(dMult-2)*(dMult-3)*(dMult-4)*(dMult-5))-9.*(pow(afvQvector1n.Mod(),4.)+pow(afvQvector2n.Mod(),2.))/(dMult*(dMult-1)*(dMult-2)*(dMult-3)*(dMult-5))+(18.*pow(afvQvector1n.Mod(),2.))/(dMult*(dMult-1)*(dMult-3)*(dMult-4))-(6.)/((dMult-1)*(dMult-2)*(dMult-3));//<6>_{n,n,n|n,n,n}
743   
744   six2n1n1n2n1n1n = (xQ2nQ1nQ1nQ2nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(2.*five2n2n2n1n1n+4.*five2n1n1n1n1n+4.*five3n1n2n1n1n+4.*four2n1n2n1n+1.*four1n1n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(4.*four1n1n1n1n+4.*two1n1n+2.*three2n1n1n+2.*three2n1n1n+4.*four3n1n1n1n+8.*three2n1n1n+2.*four4n2n1n1n+4.*four2n1n2n1n+2.*two2n2n+8.*four2n1n2n1n+4.*four3n1n3n1n+8.*three3n2n1n+4.*four3n1n2n2n+4.*four1n1n1n1n+4.*four2n1n2n1n+1.*four2n2n2n2n)-dMult*(dMult-1.)*(dMult-2.)*(2.*three2n1n1n+8.*two1n1n+4.*two1n1n+2.+4.*two1n1n+4.*three2n1n1n+2.*two2n2n+4.*three2n1n1n+8.*three3n2n1n+8.*two2n2n+4.*three4n3n1n+4.*two3n3n+4.*three3n2n1n+4.*two1n1n+8.*three2n1n1n+4.*two1n1n+4.*three3n2n1n+4.*three2n1n1n+2.*two2n2n+4.*three3n2n1n+2.*three4n2n2n)-dMult*(dMult-1.)*(4.*two1n1n+4.+4.*two1n1n+2.*two2n2n+1.+4.*two1n1n+4.*two2n2n+4.*two3n3n+   1.+2.*two2n2n+1.*two4n4n)-dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));//<6>_{2n,n,n|2n,n,n}
745  
746   six2n2n1n1n1n1n = (reQ2nQ2nQ1nstarQ1nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(five4n1n1n1n1n+8.*five2n1n1n1n1n+6.*five2n2n2n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(4.*four3n1n1n1n+6.*four4n2n1n1n+12.*three2n1n1n+12.*four1n1n1n1n+24.*four2n1n2n1n+4.*four3n1n2n2n+3.*four2n2n2n2n)-dMult*(dMult-1.)*(dMult-2.)*(6.*three2n1n1n+12.*three3n2n1n+4.*three4n3n1n+3.*three4n2n2n+8.*three2n1n1n+24.*two1n1n+12.*two2n2n+12.*three2n1n1n+8.*three3n2n1n+1.*three4n2n2n)-dMult*(dMult-1.)*(4.*two1n1n+6.*two2n2n+4.*two3n3n+1.*two4n4n+2.*two2n2n+8.*two1n1n+6.)-dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));//<6>_{2n,2n,n|n,n,n}
747    
748   six3n1n1n1n1n1n = (reQ3nQ1nQ1nstarQ1nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(five4n1n1n1n1n+4.*five2n1n1n1n1n+6.*five3n1n2n1n1n+4.*four3n1n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(4.*four3n1n1n1n+6.*four4n2n1n1n+6.*four1n1n1n1n+12.*three2n1n1n+12.*four2n1n2n1n+6.*four3n1n1n1n+12.*three3n2n1n+4.*four3n1n3n1n+3.*four3n1n2n2n)-dMult*(dMult-1.)*(dMult-2.)*(6.*three2n1n1n+12.*three3n2n1n+4.*three4n3n1n+3.*three4n2n2n+4.*two1n1n+12.*two1n1n+6.*three2n1n1n+12.*three2n1n1n+4.*three3n2n1n+12.*two2n2n+4.*three3n2n1n+4.*two3n3n+1.*three4n3n1n+6.*three3n2n1n)-dMult*(dMult-1.)*(4.*two1n1n+6.*two2n2n+4.*two3n3n+1.*two4n4n+1.*two1n1n+4.+6.*two1n1n+4.*two2n2n+1.*two3n3n)-dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));//<6>_{3n,n|n,n,n,n}
749    
750   fQCorrelations->Fill(23.,six1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
751   fQCorrelations->Fill(24.,six2n1n1n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
752   fQCorrelations->Fill(25.,six2n2n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));
753   fQCorrelations->Fill(26.,six3n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
754
755   f6pDistribution->Fill(six1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
756   
757   fQProduct->Fill(1.,two1n1n*six1n1n1n1n1n1n,dMult*(dMult-1.)*dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));
758   fQProduct->Fill(3.,four1n1n1n1n*six1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));
759  }
760  
761  //7-particle
762  Double_t seven2n1n1n1n1n1n1n=0.;
763  if(dMult>6)
764  {
765   seven2n1n1n1n1n1n1n = (reQ2nQ1nQ1nQ1nstarQ1nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(2.*six3n1n1n1n1n1n+4.*six1n1n1n1n1n1n+1.*six2n2n1n1n1n1n+6.*six2n1n1n2n1n1n+8.*five2n1n1n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(1.*five4n1n1n1n1n +8.*five2n1n1n1n1n+8.*four3n1n1n1n+12.*five3n1n2n1n1n+4.*five2n1n1n1n1n+3.*five2n2n2n1n1n+6.*five2n2n2n1n1n+6.*four1n1n1n1n+24.*four1n1n1n1n+12.*five2n1n1n1n1n+12.*five2n1n1n1n1n+12.*three2n1n1n+24.*four2n1n2n1n+4.*five3n1n2n1n1n+4.*five2n1n1n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(4.*four3n1n1n1n+6.*four4n2n1n1n+12.*four1n1n1n1n+24.*three2n1n1n+24.*four2n1n2n1n+12.*four3n1n1n1n+24.*three3n2n1n+8.*four3n1n3n1n+6.*four3n1n2n2n+6.*three2n1n1n+12.*four1n1n1n1n+12.*four2n1n2n1n+6.*three2n1n1n+12.*four2n1n2n1n+4.*four3n1n2n2n+3.*four2n2n2n2n+4.*four1n1n1n1n+6.*three2n1n1n+24.*two1n1n+24.*four1n1n1n1n+4.*four3n1n1n1n+24.*two1n1n+24.*three2n1n1n+12.*two2n2n+24.*three2n1n1n+12.*four2n1n2n1n+8.*three3n2n1n+8.*four2n1n2n1n+1.*four4n2n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(6.*three2n1n1n+1.*three2n1n1n+8.*two1n1n+12.*three3n2n1n+24.*two1n1n+12.*three2n1n1n+4.*three2n1n1n+8.*two1n1n+4.*three4n3n1n+24.*three2n1n1n+8.*three3n2n1n+12.*two1n1n+12.*two1n1n+3.*three4n2n2n+24.*two2n2n+6.*two2n2n+12.+12.*three3n2n1n+8.*two3n3n+12.*three2n1n1n+24.*two1n1n+4.*three3n2n1n+8.*three3n2n1n+2.*three4n3n1n+12.*two1n1n+8.*three2n1n1n+4.*three2n1n1n+2.*three3n2n1n+6.*two2n2n+8.*two2n2n+1.*three4n2n2n+4.*three3n2n1n+6.*three2n1n1n)-dMult*(dMult-1.)*(4.*two1n1n+2.*two1n1n+6.*two2n2n+8.+1.*two2n2n+4.*two3n3n+12.*two1n1n+4.*two1n1n+1.*two4n4n+8.*two2n2n+6.+2.*two3n3n+4.*two1n1n+1.*two2n2n)-dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.));
766         
767   fQCorrelations->Fill(28.,seven2n1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.));
768  }
769  
770  //8-particle
771  Double_t eight1n1n1n1n1n1n1n1n=0.;
772  if(dMult>7)
773  {
774   //fill the common control histogram (8th order):
775   fCommonHists8th->FillControlHistograms(anEvent); 
776   
777   eight1n1n1n1n1n1n1n1n = (pow(afvQvector1n.Mod(),8)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(12.*seven2n1n1n1n1n1n1n+16.*six1n1n1n1n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(8.*six3n1n1n1n1n1n+48.*six1n1n1n1n1n1n+6.*six2n2n1n1n1n1n+96.*five2n1n1n1n1n+72.*four1n1n1n1n+36.*six2n1n1n2n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(2.*five4n1n1n1n1n+32.*five2n1n1n1n1n+36.*four1n1n1n1n+32.*four3n1n1n1n+48.*five2n1n1n1n1n+48.*five3n1n2n1n1n+144.*five2n1n1n1n1n+288.*four1n1n1n1n+36.*five2n2n2n1n1n+144.*three2n1n1n+96.*two1n1n+144.*four2n1n2n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(8.*four3n1n1n1n+48.*four1n1n1n1n+12.*four4n2n1n1n+96.*four2n1n2n1n+96.*three2n1n1n+72.*three2n1n1n+144.*two1n1n+16.*four3n1n3n1n+48.*four3n1n1n1n+144.*four1n1n1n1n+72.*four1n1n1n1n+96.*three3n2n1n+24.*four3n1n2n2n+144.*four2n1n2n1n+288.*two1n1n+288.*three2n1n1n+9.*four2n2n2n2n+72.*two2n2n+24.)-dMult*(dMult-1.)*(dMult-2.)*(12.*three2n1n1n+16.*two1n1n+24.*three3n2n1n+48.*three2n1n1n+96.*two1n1n+8.*three4n3n1n+32.*three3n2n1n+96.*three2n1n1n+144.*two1n1n+6.*three4n2n2n+96.*two2n2n+36.*two2n2n+72.+48.*three3n2n1n+16.*two3n3n+72.*three2n1n1n+144.*two1n1n)-dMult*(dMult-1.)*(8.*two1n1n+12.*two2n2n+16.+8.*two3n3n+48.*two1n1n+1.*two4n4n+16.*two2n2n+18.)-dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.));
778   
779   fQCorrelations->Fill(30.,eight1n1n1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.));
780  
781   f8pDistribution->Fill(eight1n1n1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.));
782  } 
783  //---------------------------------------------------------------------------------------------------------
784  
785  
786  
787  
788  //---------------------------------------------------------------------------------------------------------
789  // weights:
790  Bool_t useWeights = fUsePhiWeights||fUsePtWeights||fUseEtaWeights;
791
792  TH1F *phiWeights = NULL; // histogram with phi weights
793  TH1D *ptWeights  = NULL; // histogram with pt weights
794  TH1D *etaWeights = NULL; // histogram with eta weights
795    
796  if(useWeights)
797  {
798   if(!fWeightsList)
799   {
800    cout<<" WARNING: fWeightsList is NULL pointer. "<<endl;
801    exit(0);
802   }
803   if(fUsePhiWeights) 
804   {
805    phiWeights = dynamic_cast<TH1F *>(fWeightsList->FindObject("phi_weights"));
806    if(!phiWeights)
807    {
808     cout<<" WARNING: couldn't access the histogram with phi weights. "<<endl;
809     exit(0);
810    } 
811   } 
812   if(fUsePtWeights) 
813   { 
814    ptWeights = dynamic_cast<TH1D *>(fWeightsList->FindObject("pt_weights"));
815    if(!ptWeights) 
816    {
817     cout<<" WARNING: couldn't access the histogram with pt weights. "<<endl;
818     exit(0);
819    } 
820   } 
821   if(fUseEtaWeights) 
822   {
823    etaWeights = dynamic_cast<TH1D *>(fWeightsList->FindObject("eta_weights"));
824    if(!etaWeights) 
825    {
826     cout<<" WARNING: couldn't access the histogram with eta weights. "<<endl;
827     exit(0);
828    }
829   } 
830  } 
831   
832  Int_t nBinsPhi = 0; 
833  Double_t dBinWidthPt=0.;
834  Double_t dBinWidthEta=0.;
835   
836  if(fnBinsPt)
837  {
838   dBinWidthPt=(fPtMax-fPtMin)/fnBinsPt;  
839  } 
840  if(fnBinsEta)
841  {
842   dBinWidthEta=(fEtaMax-fEtaMin)/fnBinsEta;  
843  } 
844  
845  if(fWeightsList)
846   {
847    if(fUsePhiWeights)
848    {
849     if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
850    }          
851    if(fUsePtWeights)
852    {
853     if(ptWeights)
854     {
855      Double_t dBinWidthPtW = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
856      if(dBinWidthPtW != dBinWidthPt)
857      {
858       cout<<" WARNING: dBinWidthPtW != dBinWidthPt in AFAWQC::Make()"<<endl;
859       exit(0);
860      }
861      Double_t dPtMinW = (ptWeights->GetXaxis())->GetXmin();
862      if(dPtMinW != fPtMin)
863      {
864       cout<<" WARNING: dPtMinW != fPtMin in AFAWQC::Make()"<<endl;
865       exit(0);
866      }
867     } 
868    }       
869    if(fUseEtaWeights)
870    {
871     if(etaWeights)
872     {
873      Double_t dBinWidthEtaW = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
874      if(dBinWidthEtaW != dBinWidthEta)
875      {
876       cout<<" WARNING: dBinWidthEtaW != dBinWidthEta in AFAWQC::Make()"<<endl;
877       exit(0);
878      }
879      Double_t dEtaMinW = (etaWeights->GetXaxis())->GetXmin();
880      if(dEtaMinW != fEtaMin)
881      {
882       cout<<" WARNING: dEtaMinW != fEtaMin in AFAWQC::Make()"<<endl;
883       exit(0);
884      }
885     } 
886    }          
887   } // end of if(weightsList)
888  
889  Double_t wPhi = 1.; // phi weight
890  Double_t wPt  = 1.; // pt weight
891  Double_t wEta = 1.; // eta weight
892  
893  Double_t dPhi = 0.;
894  Double_t dPt  = 0.;
895  Double_t dEta = 0.;
896
897  Double_t dQnkX[4][8] = {{0.}}; // sum_{i=1}^{M} w_i^k cos(n phi_i)
898  Double_t dQnkY[4][8] = {{0.}}; // sum_{i=1}^{M} w_i^k sin(n phi_i)
899  Double_t dSnk[4][8] = {{0.}};  //(sum_{i=1}^{M} w_i^k)^n
900  
901   for(Int_t i=0;i<nPrim;i++) // loop over all particles
902   { 
903    fTrack=anEvent->GetTrack(i);   
904    if(fTrack)
905    {
906     if(fTrack->UseForIntegratedFlow()) 
907     {
908      dPhi = fTrack->Phi();
909      dPt  = fTrack->Pt();
910      dEta = fTrack->Eta();
911      
912      // determine Phi weight: 
913      if(phiWeights && nBinsPhi)
914      {
915       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
916      }
917      // determine pt weight:    
918      if(ptWeights && dBinWidthPt)
919      {
920       wPt = ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/dBinWidthPt))); 
921      }            
922      // determine eta weight:    
923      if(etaWeights && dBinWidthEta)
924      {
925       wEta = etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/dBinWidthEta))); 
926      } 
927
928      // (Q_{n,k})_x, (Q_{n,k})_y and S_{n,k}
929      for(Int_t nn=0;nn<4;nn++)
930      {
931       for(Int_t k=0;k<8;k++)
932       {
933        dQnkX[nn][k]+=pow(wPhi*wPt*wEta,k+1)*TMath::Cos(2*(nn+1)*dPhi);
934        dQnkY[nn][k]+=pow(wPhi*wPt*wEta,k+1)*TMath::Sin(2*(nn+1)*dPhi);
935        dSnk[nn][k]+=pow(wPhi*wPt*wEta,k+1);
936       } 
937      }
938      
939     } // end of if (pTrack->UseForIntegratedFlow())
940    } // end of if (pTrack)
941    else {cerr << "no particle!!!"<<endl;}
942   } // loop over particles
943   
944   for(Int_t nn=0;nn<4;nn++)
945   {
946    for(Int_t k=0;k<8;k++)
947    {
948     dSnk[nn][k]=pow(dSnk[nn][k],nn+1);
949    }  
950   }
951   
952  //.............................................................................................. 
953  // Ms (introduced in order to simplify some Eqs. bellow)
954  Double_t dM11 = dSnk[1][0]-dSnk[0][1]; // dM11 = sum_{i,j=1,i!=j}^M w_i w_j
955  Double_t dM22 = dSnk[1][1]-dSnk[0][3]; // dM22 = sum_{i,j=1,i!=j}^M w_i^2 w_j^2
956  Double_t dM33 = dSnk[1][2]-dSnk[0][5]; // dM33 = sum_{i,j=1,i!=j}^M w_i^3 w_j^3
957  Double_t dM44 = dSnk[1][3]-dSnk[0][7]; // dM44 = sum_{i,j=1,i!=j}^M w_i^4 w_j^4
958  Double_t dM31 = dSnk[0][2]*dSnk[0][0]-dSnk[0][3]; // dM31 = sum_{i,j=1,i!=j}^M w_i^3 w_j
959  Double_t dM211 = dSnk[0][1]*dSnk[1][0]-2.*dSnk[0][2]*dSnk[0][0]-dSnk[1][1]+2.*dSnk[0][3]; // dM211 = sum_{i,j,k=1,i!=j!=k}^M w_i^2 w_j w_k
960  Double_t dM1111 = dSnk[3][0]-6.*dM211-4.*dM31-3.*dM22-dSnk[0][3]; // dM1111 = sum_{i,j,k,l=1,i!=j!=k!=l}^M w_i w_j w_k w_l
961  //..............................................................................................
962
963  //---------------------------------------------------------------------------------------------------------
964  //
965  //                                        ***********************************************
966  //                                        **** weighted multi-particle correlations: ****
967  //                                        ***********************************************
968  //
969  // Remark 1: weighted multi-particle correlations calculated with Q-vectors are stored in fWeightedQCorrelations.
970  // Remark 2: binning of fWeightedQCorrelations is organized as follows: 
971  
972  // binning
973  //..............................................................................................
974  // 1st bin: weighted <2>_{n|n}   = <w1   w2   cos( n*(phi1-phi2))>
975  // 2nd bin: weighted <2>_{2n|2n} = <w1^2 w2^2 cos(2n*(phi1-phi2))>
976  // 3rd bin: weighted <2>_{3n|3n} = <w1^3 w2^3 cos(3n*(phi1-phi2))>
977  // 4th bin: weighted <2>_{4n|4n} = <w1^4 w2^4 cos(4n*(phi1-phi2))>
978  // 5th bin: weighted <2>_{n|n} = <w1^3 w2 cos(n*(phi1-phi2))>
979  // 6th bin: weighted <2>_{n|n} = <w1 w2 w3^2 cos(n*(phi1-phi2))>
980  
981  // 11th bin: weighted <3>_{2n|n,n} = <w1^2 w2 w3 cos(n*(2phi1-phi2-phi3))> 
982  
983  // 21st bin: weighted <4>_{n,n|n,n} = <w1 w2 w3 w4 cos(n*(phi1+phi2-phi3-phi4))>
984  //..............................................................................................
985  
986  //.............................................................................................. 
987  // weighted 2-particle correlations:
988  Double_t two1n1nW1W1=0., two2n2nW2W2=0., two3n3nW3W3=0., two4n4nW4W4=0., two1n1nW3W1=0., two1n1nW2W1W1=0.;
989  
990  if(nRP>1) // nRP = number of particles used to determine the reaction plane
991  { 
992   if(dM11 != 0)
993   {
994    two1n1nW1W1 = (pow(dQnkX[0][0],2)+pow(dQnkY[0][0],2)-dSnk[0][1])/dM11; // <2>_{n|n}=<w1 w2 cos(n*(phi1-phi2))>
995    fWeightedQCorrelations->Fill(0.,two1n1nW1W1,dM11);
996   }
997   if(dM22 != 0)
998   {
999    two2n2nW2W2 = (pow(dQnkX[1][1],2)+pow(dQnkY[1][1],2)-dSnk[0][3])/dM22; // <2>_{2n|2n}=<w1^2 w2^2 cos(2n*(phi1-phi2))>
1000    fWeightedQCorrelations->Fill(1.,two2n2nW2W2,dM22); 
1001   }
1002   if(dM33 != 0)
1003   {
1004    two3n3nW3W3 = (pow(dQnkX[2][2],2)+pow(dQnkY[2][2],2)-dSnk[0][5])/dM33; // <2>_{2n|2n}=<w1^3 w2^3 cos(3n*(phi1-phi2))>
1005    fWeightedQCorrelations->Fill(2.,two3n3nW3W3,dM33);
1006   }
1007   if(dM44 != 0)
1008   {
1009    two4n4nW4W4 = (pow(dQnkX[3][3],2)+pow(dQnkY[3][3],2)-dSnk[0][7])/dM44; // <2>_{2n|2n}=<w1^4 w2^4 cos(4n*(phi1-phi2))>
1010    fWeightedQCorrelations->Fill(3.,two4n4nW4W4,dM44);  
1011   } 
1012   if(dM31 != 0)
1013   {
1014    two1n1nW3W1 = (dQnkX[0][2]*dQnkX[0][0]+dQnkY[0][2]*dQnkY[0][0]-dSnk[0][3])/dM31; // <2>_{n|n}=<w1^3 w2 cos(n*(phi1-phi2))>
1015    fWeightedQCorrelations->Fill(4.,two1n1nW3W1,dM31);  
1016   } 
1017   if(dM211 != 0)
1018   {
1019    two1n1nW2W1W1 = (dSnk[0][1]*dM11*two1n1nW1W1-2.*dM31*two1n1nW3W1)/dM211; // <2>_{n|n}=<w1^2 w2 w3 cos(n*(phi1-phi2))>
1020    fWeightedQCorrelations->Fill(5.,two1n1nW2W1W1,dM211);  
1021   } 
1022  } // end of if(nRP>1)
1023  //..............................................................................................
1024
1025  //..............................................................................................
1026  // weighted 3-particle correlations:
1027  Double_t three2n1n1nW2W1W1=0.;
1028  
1029  if(nRP>2) // nRP = number of particles used to determine the reaction plane
1030  { 
1031   if(dM211 != 0)
1032   {
1033    three2n1n1nW2W1W1 = (pow(dQnkX[0][0],2.)*dQnkX[1][1]+2.*dQnkX[0][0]*dQnkY[0][0]*dQnkY[1][1]-pow(dQnkY[0][0],2.)*dQnkX[1][1]-2.*dM31*two1n1nW3W1-dM22*two2n2nW2W2-dSnk[0][3])/dM211;
1034    fWeightedQCorrelations->Fill(10.,three2n1n1nW2W1W1,dM211);
1035   } 
1036  } // end of if(nRP>2) 
1037  //..............................................................................................
1038  
1039  //..............................................................................................
1040  // weighted 4-particle correlations:
1041  Double_t four1n1n1n1nW1W1W1W1=0.;
1042  
1043  if(nRP>3) // nRP = number of particles used to determine the reaction plane
1044  { 
1045   if(dM1111 != 0)
1046   {
1047    four1n1n1n1nW1W1W1W1 = (pow(pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.),2)-2.*dM211*three2n1n1nW2W1W1-4.*dM211*two1n1nW2W1W1-4.*dM31*two1n1nW3W1-dM22*two2n2nW2W2-2.*dM22-dSnk[0][3])/dM1111;
1048    fWeightedQCorrelations->Fill(20.,four1n1n1n1nW1W1W1W1,dM1111);
1049   } 
1050  } // end of if(nRP>3) 
1051  //..............................................................................................
1052      
1053  
1054  
1055  
1056  //---------------------------------------------------------------------------------------------------------
1057  //
1058  //                                        *******************************************************
1059  //                                        **** weighted reduced multi-particle correlations: ****
1060  //                                        *******************************************************
1061  // 
1062  // pt POI
1063  TProfile *ptReq1nPrime = new TProfile("ptReq1nPrime","Re[q_{n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1064  TProfile *ptImq1nPrime = new TProfile("ptImq1nPrime","Im[q_{n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1065  TProfile *ptReq2nPrime = new TProfile("ptReq2nPrime","Re[q_{2n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1066  TProfile *ptImq2nPrime = new TProfile("ptImq2nPrime","Im[q_{2n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1067  
1068  TProfile *ptReq1nPrimePrime = new TProfile("ptReq1nPrimePrime","Re[q_{n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1069  TProfile *ptImq1nPrimePrime = new TProfile("ptImq1nPrimePrime","Im[q_{n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1070  TProfile *ptReq2nPrimePrime = new TProfile("ptReq2nPrimePrime","Re[q_{2n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1071  TProfile *ptImq2nPrimePrime = new TProfile("ptImq2nPrimePrime","Im[q_{2n}^{''}]",fnBinsPt,fPtMin,fPtMax,"s");
1072  
1073  TProfile *req1nW2PrimePrimePt = new TProfile("req1nW2PrimePrimePt","#sum_{i=1}^{m''} w_{i}^{2} cos(n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1074  TProfile *imq1nW2PrimePrimePt = new TProfile("imq1nW2PrimePrimePt","#sum_{i=1}^{m''} w_{i}^{2} sin(n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1075  TProfile *req2nW1PrimePrimePt = new TProfile("req2nW1PrimePrimePt","#sum_{i=1}^{m''} w_{i} cos(2n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1076  TProfile *imq2nW1PrimePrimePt = new TProfile("imq2nW1PrimePrimePt","#sum_{i=1}^{m''} w_{i} sin(2n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1077  
1078  TProfile *sumOfW1upTomPrimePrimePt = new TProfile("sumOfW1upTomPrimePrimePt","#sum_{i=1}^{m''} w_{i}''",fnBinsPt,fPtMin,fPtMax,"s");
1079  TProfile *sumOfW2upTomPrimePrimePt = new TProfile("sumOfW2upTomPrimePrimePt","#sum_{i=1}^{m''} w_{i}^{2}''",fnBinsPt,fPtMin,fPtMax,"s");
1080  TProfile *sumOfW3upTomPrimePrimePt = new TProfile("sumOfW3upTomPrimePrimePt","#sum_{i=1}^{m''} w_{i}^{3}''",fnBinsPt,fPtMin,fPtMax,"s");
1081  
1082  // eta POI 
1083  TProfile *etaReq1nPrime = new TProfile("etaReq1nPrime","Re[q_{n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1084  TProfile *etaImq1nPrime = new TProfile("etaImq1nPrime","Im[q_{n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1085  TProfile *etaReq2nPrime = new TProfile("etaReq2nPrime","Re[q_{2n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1086  TProfile *etaImq2nPrime = new TProfile("etaImq2nPrime","Im[q_{2n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1087  
1088  TProfile *etaReq1nPrimePrime = new TProfile("etaReq1nPrimePrime","Re[q_{n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1089  TProfile *etaImq1nPrimePrime = new TProfile("etaImq1nPrimePrime","Im[q_{n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1090  TProfile *etaReq2nPrimePrime = new TProfile("etaReq2nPrimePrime","Re[q_{2n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1091  TProfile *etaImq2nPrimePrime = new TProfile("etaImq2nPrimePrime","Im[q_{2n}^{''}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1092  
1093  TProfile *req1nW2PrimePrimeEta = new TProfile("req1nW2PrimePrimeEta","#sum_{i=1}^{m''} w_{i}^{2} cos(n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1094  TProfile *imq1nW2PrimePrimeEta = new TProfile("imq1nW2PrimePrimeEta","#sum_{i=1}^{m''} w_{i}^{2} sin(n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1095  TProfile *req2nW1PrimePrimeEta = new TProfile("req2nW1PrimePrimeEta","#sum_{i=1}^{m''} w_{i} cos(2n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1096  TProfile *imq2nW1PrimePrimeEta = new TProfile("imq2nW1PrimePrimeEta","#sum_{i=1}^{m''} w_{i} sin(2n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1097  
1098  TProfile *sumOfW1upTomPrimePrimeEta = new TProfile("sumOfW1upTomPrimePrimeEta","#sum_{i=1}^{m''} w_{i}''",fnBinsEta,fEtaMin,fEtaMax,"s");
1099  TProfile *sumOfW2upTomPrimePrimeEta = new TProfile("sumOfW2upTomPrimePrimeEta","#sum_{i=1}^{m''} w_{i}^{2}''",fnBinsEta,fEtaMin,fEtaMax,"s");
1100  TProfile *sumOfW3upTomPrimePrimeEta = new TProfile("sumOfW3upTomPrimePrimeEta","#sum_{i=1}^{m''} w_{i}^{3}''",fnBinsEta,fEtaMin,fEtaMax,"s");
1101  
1102  // pt RP
1103  TProfile *ptReq1n = new TProfile("ptReq1n","Re[q_{n}]",fnBinsPt,fPtMin,fPtMax,"s");
1104  TProfile *ptImq1n = new TProfile("ptImq1n","Im[q_{n}]",fnBinsPt,fPtMin,fPtMax,"s");
1105  TProfile *ptReq2n = new TProfile("ptReq2n","Re[q_{2n}]",fnBinsPt,fPtMin,fPtMax,"s");
1106  TProfile *ptImq2n = new TProfile("ptImq2n","Im[q_{2n}]",fnBinsPt,fPtMin,fPtMax,"s");
1107  
1108  TProfile *req1nW2Pt = new TProfile("req1nW2Pt","#sum_{i=1}^{m} w_{i}^{2} cos(n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1109  TProfile *imq1nW2Pt = new TProfile("imq1nW2Pt","#sum_{i=1}^{m} w_{i}^{2} sin(n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1110  TProfile *req2nW1Pt = new TProfile("req2nW1Pt","#sum_{i=1}^{m} w_{i} cos(2n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1111  TProfile *imq2nW1Pt = new TProfile("imq2nW1Pt","#sum_{i=1}^{m} w_{i} sin(2n(#phi_{i}))''",fnBinsPt,fPtMin,fPtMax,"s");
1112  
1113  TProfile *sumOfW1upTomPt = new TProfile("sumOfW1upTomPt","#sum_{i=1}^{m} w_{i}''",fnBinsPt,fPtMin,fPtMax,"s");
1114  TProfile *sumOfW2upTomPt = new TProfile("sumOfW2upTomPt","#sum_{i=1}^{m} w_{i}^{2}''",fnBinsPt,fPtMin,fPtMax,"s");
1115  TProfile *sumOfW3upTomPt = new TProfile("sumOfW3upTomPt","#sum_{i=1}^{m} w_{i}^{3}''",fnBinsPt,fPtMin,fPtMax,"s");
1116  
1117  // eta RP
1118  TProfile *etaReq1n = new TProfile("etaReq1n","Re[q_{n}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1119  TProfile *etaImq1n = new TProfile("etaImq1n","Im[q_{n}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1120  TProfile *etaReq2n = new TProfile("etaReq2n","Re[q_{2n}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1121  TProfile *etaImq2n = new TProfile("etaImq2n","Im[q_{2n}]",fnBinsEta,fEtaMin,fEtaMax,"s");
1122  
1123  TProfile *req1nW2Eta = new TProfile("req1nW2Eta","#sum_{i=1}^{m} w_{i}^{2} cos(n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1124  TProfile *imq1nW2Eta = new TProfile("imq1nW2Eta","#sum_{i=1}^{m} w_{i}^{2} sin(n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1125  TProfile *req2nW1Eta = new TProfile("req2nW1Eta","#sum_{i=1}^{m} w_{i} cos(2n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1126  TProfile *imq2nW1Eta = new TProfile("imq2nW1Eta","#sum_{i=1}^{m} w_{i} sin(2n(#phi_{i}))''",fnBinsEta,fEtaMin,fEtaMax,"s");
1127  
1128  TProfile *sumOfW1upTomEta = new TProfile("sumOfW1upTomEta","#sum_{i=1}^{m} w_{i}''",fnBinsEta,fEtaMin,fEtaMax,"s");
1129  TProfile *sumOfW2upTomEta = new TProfile("sumOfW2upTomEta","#sum_{i=1}^{m} w_{i}^{2}''",fnBinsEta,fEtaMin,fEtaMax,"s");
1130  TProfile *sumOfW3upTomEta = new TProfile("sumOfW3upTomEta","#sum_{i=1}^{m} w_{i}^{3}''",fnBinsEta,fEtaMin,fEtaMax,"s");
1131
1132  for(Int_t i=0;i<nPrim;i++) // loop over all particles  
1133  { 
1134   fTrack=anEvent->GetTrack(i);
1135   if(fTrack)
1136   {
1137    if(fTrack->UseForDifferentialFlow()) // checking if particle is POI 
1138    {
1139     if(fTrack->UseForIntegratedFlow()) // checking if particle is both POI and RP 
1140     {
1141      // get azimuthal angle, momentum and pseudorapidity of a particle:
1142      dPhi = fTrack->Phi();
1143      dPt  = fTrack->Pt();
1144      dEta = fTrack->Eta();
1145      // phi weights:
1146      if(fUsePhiWeights) 
1147      {
1148       nBinsPhi = phiWeights->GetNbinsX();
1149       if(nBinsPhi) 
1150       {
1151        wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
1152       }
1153      } 
1154      // pt weights:
1155      if(fUsePtWeights)
1156      {          
1157       if(dBinWidthPt) 
1158       {
1159        wPt = ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/dBinWidthPt)));
1160       }
1161      }             
1162      // eta weights:
1163      if(fUseEtaWeights)
1164      {    
1165       if(dBinWidthEta)
1166       {
1167        wEta = etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/dBinWidthEta))); 
1168       }
1169      }
1170      // pt:
1171      ptReq1nPrimePrime->Fill(dPt,cos(n*dPhi),1.); 
1172      ptImq1nPrimePrime->Fill(dPt,sin(n*dPhi),1.);
1173      ptReq2nPrimePrime->Fill(dPt,cos(2.*n*dPhi),1.);
1174      ptImq2nPrimePrime->Fill(dPt,sin(2.*n*dPhi),1.);
1175      // weighted pt:
1176      req1nW2PrimePrimePt->Fill(dPt,cos(n*dPhi),pow(wPhi*wPt*wEta,2.));
1177      imq1nW2PrimePrimePt->Fill(dPt,sin(n*dPhi),pow(wPhi*wPt*wEta,2.));
1178      req2nW1PrimePrimePt->Fill(dPt,cos(2.*n*dPhi),wPhi*wPt*wEta);
1179      imq2nW1PrimePrimePt->Fill(dPt,sin(2.*n*dPhi),wPhi*wPt*wEta);
1180      sumOfW1upTomPrimePrimePt->Fill(dPt,wPhi*wPt*wEta,1.);
1181      sumOfW2upTomPrimePrimePt->Fill(dPt,pow(wPhi*wPt*wEta,2),1.);
1182      sumOfW3upTomPrimePrimePt->Fill(dPt,pow(wPhi*wPt*wEta,3),1.);
1183      
1184      // eta:
1185      etaReq1nPrimePrime->Fill(dEta,cos(n*dPhi),1.); 
1186      etaImq1nPrimePrime->Fill(dEta,sin(n*dPhi),1.);
1187      etaReq2nPrimePrime->Fill(dEta,cos(2.*n*dPhi),1.);
1188      etaImq2nPrimePrime->Fill(dEta,sin(2.*n*dPhi),1.);
1189      // weighted eta:
1190      req1nW2PrimePrimeEta->Fill(dEta,cos(n*dPhi),pow(wPhi*wPt*wEta,2.));
1191      imq1nW2PrimePrimeEta->Fill(dEta,sin(n*dPhi),pow(wPhi*wPt*wEta,2.));
1192      req2nW1PrimePrimeEta->Fill(dEta,cos(2.*n*dPhi),wPhi*wPt*wEta);
1193      imq2nW1PrimePrimeEta->Fill(dEta,sin(2.*n*dPhi),wPhi*wPt*wEta);
1194      sumOfW1upTomPrimePrimeEta->Fill(dEta,wPhi*wPt*wEta,1.);
1195      sumOfW2upTomPrimePrimeEta->Fill(dEta,pow(wPhi*wPt*wEta,2),1.);
1196      sumOfW3upTomPrimePrimeEta->Fill(dEta,pow(wPhi*wPt*wEta,3),1.);
1197
1198     }else if(!(fTrack->UseForIntegratedFlow())) // checking if particles is POI and not RP  
1199      {
1200       // get azimuthal angle, momentum and pseudorapidity of a particle:
1201       dPhi = fTrack->Phi();
1202       dPt  = fTrack->Pt();
1203       dEta = fTrack->Eta();
1204       // pt:
1205       ptReq1nPrime->Fill(dPt,cos(n*dPhi),1.);   
1206       ptImq1nPrime->Fill(dPt,sin(n*dPhi),1.);
1207       ptReq2nPrime->Fill(dPt,cos(2.*n*dPhi),1.);
1208       ptImq2nPrime->Fill(dPt,sin(2.*n*dPhi),1.);
1209
1210       // eta:
1211       etaReq1nPrime->Fill(dEta,cos(n*dPhi),1.); 
1212       etaImq1nPrime->Fill(dEta,sin(n*dPhi),1.);
1213       etaReq2nPrime->Fill(dEta,cos(2.*n*dPhi),1.);
1214       etaImq2nPrime->Fill(dEta,sin(2.*n*dPhi),1.);
1215
1216      } // end of else if(!(fTrack->UseForIntegratedFlow())) // checking if particles is POI and not RP 
1217    } // end of if(fTrack->UseForDifferentialFlow()) // checking if particle is POI 
1218      
1219    if(fTrack->UseForIntegratedFlow()) // checking if particles is only RP:
1220    {
1221     dPhi = fTrack->Phi();
1222     dPt  = fTrack->Pt();
1223     dEta = fTrack->Eta();
1224      
1225     // phi weights:
1226     if(fUsePhiWeights) 
1227     {
1228      nBinsPhi = phiWeights->GetNbinsX();
1229      if(nBinsPhi) 
1230      {
1231       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
1232      }
1233     } 
1234     // pt weights:
1235     if(fUsePtWeights)
1236     {          
1237      if(dBinWidthPt) 
1238      {
1239       wPt = ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/dBinWidthPt)));
1240      }
1241     }             
1242     // eta weights:
1243     if(fUseEtaWeights)
1244     {    
1245      if(dBinWidthEta)
1246      {
1247       wEta = etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/dBinWidthEta))); 
1248      }
1249     }
1250     // pt:
1251     ptReq1n->Fill(dPt,cos(n*dPhi),1.); 
1252     ptImq1n->Fill(dPt,sin(n*dPhi),1.);
1253     ptReq2n->Fill(dPt,cos(2.*n*dPhi),1.);
1254     ptImq2n->Fill(dPt,sin(2.*n*dPhi),1.);
1255     // weighted pt:
1256     req1nW2Pt->Fill(dPt,cos(n*dPhi),pow(wPhi*wPt*wEta,2.));
1257     imq1nW2Pt->Fill(dPt,sin(n*dPhi),pow(wPhi*wPt*wEta,2.));
1258     req2nW1Pt->Fill(dPt,cos(2.*n*dPhi),wPhi*wPt*wEta);
1259     imq2nW1Pt->Fill(dPt,sin(2.*n*dPhi),wPhi*wPt*wEta);
1260     sumOfW1upTomPt->Fill(dPt,wPhi*wPt*wEta,1.);
1261     sumOfW2upTomPt->Fill(dPt,pow(wPhi*wPt*wEta,2),1.);
1262     sumOfW3upTomPt->Fill(dPt,pow(wPhi*wPt*wEta,3),1.);
1263      
1264     // eta:
1265     etaReq1n->Fill(dEta,cos(n*dPhi),1.); 
1266     etaImq1n->Fill(dEta,sin(n*dPhi),1.);
1267     etaReq2n->Fill(dEta,cos(2.*n*dPhi),1.);
1268     etaImq2n->Fill(dEta,sin(2.*n*dPhi),1.);
1269     // weighted eta:
1270     req1nW2Eta->Fill(dEta,cos(n*dPhi),pow(wPhi*wPt*wEta,2.));
1271     imq1nW2Eta->Fill(dEta,sin(n*dPhi),pow(wPhi*wPt*wEta,2.));
1272     req2nW1Eta->Fill(dEta,cos(2.*n*dPhi),wPhi*wPt*wEta);
1273     imq2nW1Eta->Fill(dEta,sin(2.*n*dPhi),wPhi*wPt*wEta);
1274     sumOfW1upTomEta->Fill(dEta,wPhi*wPt*wEta,1.);
1275     sumOfW2upTomEta->Fill(dEta,pow(wPhi*wPt*wEta,2),1.);
1276     sumOfW3upTomEta->Fill(dEta,pow(wPhi*wPt*wEta,3),1.);
1277    } // end of if(fTrack->UseForIntegratedFlow()) // checking if particles is only RP:
1278   } // end of if(fTrack}      
1279  } // end of for(Int_t i=0;i<nPrim;i++) 
1280  
1281  //...........................................................................................................
1282  // PrimePrime Pt POI
1283  Double_t qxPrimePrimePtPOI=0.,qyPrimePrimePtPOI=0.,q2xPrimePrimePtPOIHere=0.,q2yPrimePrimePtPOIHere=0.;//add comments for these variable
1284  Double_t qxW2PrimePrimePtPOI=0.,qyW2PrimePrimePtPOI=0.,q2xW1PrimePrimePtPOI=0.,q2yW1PrimePrimePtPOI=0.;//add comments for these variable
1285  Double_t dS11mPrimePrimePtPOI=0.; // to be improved (name)
1286  Double_t dS12mPrimePrimePtPOI=0.; // to be improved (name)
1287  Double_t dS13mPrimePrimePtPOI=0.; // to be improved (name)
1288  Double_t mPrimePrimePtPOI=0.; // to be improved (name)
1289
1290  Double_t dM1pp11PtPOI=0.; // to be improved (name)
1291  Double_t dM0pp111PtPOI=0.; // to be improved (name)
1292  Double_t dM0pp12PtPOI=0.;
1293  Double_t dM2pp1PtPOI=0.;
1294  Double_t dM1pp2PtPOI=0.;
1295  Double_t dM0pp3PtPOI=0.;
1296  
1297  // Prime Pt POI
1298  Double_t qxPrimePtPOI=0.,qyPrimePtPOI=0.;
1299  Double_t mPrimePtPOI=0.;
1300  Double_t dM0p111PtPOI=0.; // to be improved (name)
1301   
1302  for(Int_t bin=1;bin<(fnBinsPt+1);bin++) // loop over pt-bins 
1303  {     
1304   // q'':                    
1305   qxPrimePrimePtPOI = (ptReq1nPrimePrime->GetBinContent(bin))*(ptReq1nPrimePrime->GetBinEntries(bin));
1306   qyPrimePrimePtPOI = (ptImq1nPrimePrime->GetBinContent(bin))*(ptImq1nPrimePrime->GetBinEntries(bin)); 
1307   q2xPrimePrimePtPOIHere = (ptReq2nPrimePrime->GetBinContent(bin))*(ptReq2nPrimePrime->GetBinEntries(bin));  
1308   q2yPrimePrimePtPOIHere = (ptImq2nPrimePrime->GetBinContent(bin))*(ptImq2nPrimePrime->GetBinEntries(bin)); 
1309   
1310   qxW2PrimePrimePtPOI = (req1nW2PrimePrimePt->GetBinContent(bin))*(req1nW2PrimePrimePt->GetBinEntries(bin));
1311   qyW2PrimePrimePtPOI = (imq1nW2PrimePrimePt->GetBinContent(bin))*(imq1nW2PrimePrimePt->GetBinEntries(bin));
1312   
1313   q2xW1PrimePrimePtPOI = (req2nW1PrimePrimePt->GetBinContent(bin))*(req2nW1PrimePrimePt->GetBinEntries(bin));
1314   q2yW1PrimePrimePtPOI = (imq2nW1PrimePrimePt->GetBinContent(bin))*(imq2nW1PrimePrimePt->GetBinEntries(bin));
1315   
1316   dS11mPrimePrimePtPOI = (sumOfW1upTomPrimePrimePt->GetBinContent(bin))*(sumOfW1upTomPrimePrimePt->GetBinEntries(bin));
1317   dS12mPrimePrimePtPOI = (sumOfW2upTomPrimePrimePt->GetBinContent(bin))*(sumOfW2upTomPrimePrimePt->GetBinEntries(bin));
1318   dS13mPrimePrimePtPOI = (sumOfW3upTomPrimePrimePt->GetBinContent(bin))*(sumOfW3upTomPrimePrimePt->GetBinEntries(bin));
1319
1320   mPrimePrimePtPOI = sumOfW1upTomPrimePrimePt->GetBinEntries(bin); // to be improved
1321       
1322   dM1pp11PtPOI=dS11mPrimePrimePtPOI*(dSnk[1][0]-dSnk[0][1])-2.*dS12mPrimePrimePtPOI*dSnk[0][0]+2.*dS13mPrimePrimePtPOI;
1323   dM1pp2PtPOI=dS11mPrimePrimePtPOI*dSnk[0][1]-dS13mPrimePrimePtPOI;
1324   dM2pp1PtPOI=dS12mPrimePrimePtPOI*dSnk[0][0]-dS13mPrimePrimePtPOI;
1325   dM0pp3PtPOI=mPrimePrimePtPOI*dSnk[0][2]-dS13mPrimePrimePtPOI;
1326   dM0pp12PtPOI=mPrimePrimePtPOI*dSnk[0][0]*dSnk[0][1]-dM2pp1PtPOI-dM1pp2PtPOI-dM0pp3PtPOI-dS13mPrimePrimePtPOI;
1327   dM0pp111PtPOI=mPrimePrimePtPOI*dSnk[2][0]-3.*dM1pp11PtPOI-3.*dM0pp12PtPOI-3.*dM2pp1PtPOI-3.*dM1pp2PtPOI-dM0pp3PtPOI-dS13mPrimePrimePtPOI;
1328   
1329   // q':
1330   qxPrimePtPOI = (ptReq1nPrime->GetBinContent(bin))*(ptReq1nPrime->GetBinEntries(bin));
1331   qyPrimePtPOI = (ptImq1nPrime->GetBinContent(bin))*(ptImq1nPrime->GetBinEntries(bin));
1332   
1333   mPrimePtPOI = ptReq1nPrime->GetBinEntries(bin); // to be improved
1334   dM0p111PtPOI=mPrimePtPOI*(dSnk[2][0]-3.*dSnk[0][1]*dSnk[0][0]+2.*dSnk[0][2]);
1335  
1336   // 2-p the needed one
1337   Double_t two1n1nWPerPtBinPOI=0.; 
1338   if((mPrimePrimePtPOI+mPrimePtPOI)*dSnk[0][0]-dS11mPrimePrimePtPOI>0)
1339   {
1340    two1n1nWPerPtBinPOI = (qxPrimePrimePtPOI*dQnkX[0][0]+qyPrimePrimePtPOI*dQnkY[0][0]+qxPrimePtPOI*dQnkX[0][0]+qyPrimePtPOI*dQnkY[0][0]-dS11mPrimePrimePtPOI)/((mPrimePrimePtPOI+mPrimePtPOI)*dSnk[0][0]-dS11mPrimePrimePtPOI); 
1341   
1342    f2WPerPtBin1n1nPOI->Fill(fPtMin+(bin-1)*dBinWidthPt,two1n1nWPerPtBinPOI,(mPrimePrimePtPOI+mPrimePtPOI)*dSnk[0][0]-dS11mPrimePrimePtPOI);     
1343   }
1344  
1345   // 2-p temporary one
1346   Double_t two1n1nW1ppW1W1PtPOI=0.; // <w1 w2 w3 cos(n(phi2-phi3))> // OK!!!
1347   if(dM1pp11PtPOI)
1348   {
1349    two1n1nW1ppW1W1PtPOI = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*dS11mPrimePrimePtPOI-2.*(qxW2PrimePrimePtPOI*dQnkX[0][0]+qyW2PrimePrimePtPOI*dQnkY[0][0]) - dS11mPrimePrimePtPOI*dSnk[0][1]+2.*dS13mPrimePrimePtPOI)/dM1pp11PtPOI; // CORRECT !!! <w1 w2 w3 cos(n(phi2-phi3))>
1350   }
1351   
1352   // 2-p temporary one
1353   Double_t two1npp1nW1W2PtPOI=0.; // <w2 w3^2 cos(n(psi1-phi2))> // OK !!!
1354   if(dM0pp12PtPOI)
1355   {
1356    two1npp1nW1W2PtPOI = (dSnk[0][1]*(qxPrimePrimePtPOI*dQnkX[0][0]+qyPrimePrimePtPOI*dQnkY[0][0])-(qxW2PrimePrimePtPOI*dQnkX[0][0]+qyW2PrimePrimePtPOI*dQnkY[0][0])-dM1pp2PtPOI-(qxPrimePrimePtPOI*dQnkX[0][2]+qyPrimePrimePtPOI*dQnkY[0][2])+dS13mPrimePrimePtPOI)/dM0pp12PtPOI; // CORRECT !!! <w2 w3^2 cos(n(psi1-phi2))>
1357   }
1358   
1359   // 2-p temporary one
1360   Double_t two1npp1nW2ppW1PtPOI=0.; // <w1^2 w2 cos(n(psi1-phi2))> // OK !!!
1361   if(dM2pp1PtPOI)
1362   {
1363    two1npp1nW2ppW1PtPOI = ((qxW2PrimePrimePtPOI*dQnkX[0][0]+qyW2PrimePrimePtPOI*dQnkY[0][0])-dS13mPrimePrimePtPOI)/dM2pp1PtPOI; // CORRECT !!! <w1^2 w2 cos(n(psi1-phi2))> 
1364   }
1365   
1366   // 2-p temporary one
1367   Double_t two2npp2nW1ppW2PtPOI=0.; // <w1 w2^2 cos(2n(psi1-phi2))> OK !!!
1368   if(dM1pp2PtPOI)
1369   {
1370    two2npp2nW1ppW2PtPOI = ((q2xW1PrimePrimePtPOI*dQnkX[1][1]+q2yW1PrimePrimePtPOI*dQnkY[1][1])-dS13mPrimePrimePtPOI)/dM1pp2PtPOI; // CORRECT !!! <w1 w2^2 cos(2n(psi1-phi2))> 
1371   }
1372   
1373   // 2-p temporary one
1374   Double_t two1npp1nW3PtPOI=0.; // <w2^3 cos(n(psi1-phi2))> // OK !!!
1375   if(dM0pp3PtPOI)
1376   {
1377    two1npp1nW3PtPOI = (qxPrimePrimePtPOI*dQnkX[0][2]+qyPrimePrimePtPOI*dQnkY[0][2]-dS13mPrimePrimePtPOI)/dM0pp3PtPOI; // CORRECT !!! <w2^3 cos(n(psi1-phi2))>
1378   }
1379    
1380   // 3-p temporary one
1381   Double_t three2npp1n1nW1ppW1W1PtPOI=0.; // <w1 w2 w3 cos(n(2psi1-phi2-phi3))> // OK!!!
1382   if(dM1pp11PtPOI)
1383   {
1384    three2npp1n1nW1ppW1W1PtPOI = (q2xW1PrimePrimePtPOI*(dQnkX[0][0]*dQnkX[0][0]-dQnkY[0][0]*dQnkY[0][0])+2.*q2yW1PrimePrimePtPOI*dQnkX[0][0]*dQnkY[0][0]-2.*(qxW2PrimePrimePtPOI*dQnkX[0][0]+qyW2PrimePrimePtPOI*dQnkY[0][0])-(q2xW1PrimePrimePtPOI*dQnkX[1][1]+q2yW1PrimePrimePtPOI*dQnkY[1][1])+2.*dS13mPrimePrimePtPOI)/dM1pp11PtPOI; // CORRECT !!! <w1 w2 w3 cos(n(2psi1-phi2-phi3))> 
1385   }
1386   
1387   // 3-p temporary one
1388   Double_t three1npp1n2nW0ppW1W2PtPOI=0.; // <w2 w3^2 cos(n(psi1+phi2-2*phi3))> // OK!!!
1389   if(dM0pp12PtPOI)
1390   {
1391    three1npp1n2nW0ppW1W2PtPOI = (qxPrimePrimePtPOI*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])-qyPrimePrimePtPOI*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1])-(qxW2PrimePrimePtPOI*dQnkX[0][0]+qyW2PrimePrimePtPOI*dQnkY[0][0])-(q2xW1PrimePrimePtPOI*dQnkX[1][1]+q2yW1PrimePrimePtPOI*dQnkY[1][1])-(qxPrimePrimePtPOI*dQnkX[0][2]+qyPrimePrimePtPOI*dQnkY[0][2])+2.*dS13mPrimePrimePtPOI)/dM0pp12PtPOI; // CORRECT !!! <w2 w3^2 cos(n(psi1+phi2-2.*phi3))>
1392   }
1393  
1394   /*
1395   // 4-p RP part
1396   Double_t four1npp1n1n1nW1W1W1=0.; // <w1 w2 w3 cos(n(psi1+phi1-phi2-phi3))>
1397   if(dM0pp111PtPOI)
1398   {   
1399    four1npp1n1n1nW1W1W1 = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePrimePtPOI*dQnkX[0][0]+qyPrimePrimePtPOI*dQnkY[0][0])-2.*dM1pp11PtPOI*two1n1nW1ppW1W1PtPOI-dM1pp11PtPOI*three2npp1n1nW1ppW1W1PtPOI-dM0pp12PtPOI*three1npp1n2nW0ppW1W2PtPOI-2.*dM0pp12PtPOI*two1npp1nW1W2PtPOI-3.*dM2pp1PtPOI*two1npp1nW2ppW1PtPOI-2.*dM1pp2PtPOI-dM1pp2PtPOI*two2npp2nW1ppW2PtPOI-dM0pp3PtPOI*two1npp1nW3PtPOI-dS13mPrimePrimePtPOI)/(dM0pp111PtPOI); 
1400   }
1401   */
1402   
1403   /*
1404   // 4-p POI part
1405   Double_t four1npp1n1n1nW1W1W1POI=0.;
1406   if(dM0p111PtPOI>0&&mPrimePtPOI>0&&nRP>0)
1407   {
1408    four1npp1n1n1nW1W1W1POI = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePtPOI*dQnkX[0][0]+qyPrimePtPOI*dQnkY[0][0])-2.*dSnk[0][1]* (qxPrimePtPOI*dQnkX[0][0]+qyPrimePtPOI*dQnkY[0][0])+2.*(qxPrimePtPOI*dQnkX[0][2]+qyPrimePtPOI*dQnkY[0][2])-qxPrimePtPOI*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimePtPOI*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/dM0p111PtPOI;
1409   }
1410   */
1411  
1412   // 4-p RP and POI in all combinations (full, partial and no overlap)
1413   Double_t four1npp1n1n1nW1W1W1PtPOI=0.;
1414   if(dM0pp111PtPOI+dM0p111PtPOI)
1415   {
1416   four1npp1n1n1nW1W1W1PtPOI = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePrimePtPOI*dQnkX[0][0]+qyPrimePrimePtPOI*dQnkY[0][0])-2.*dM1pp11PtPOI*two1n1nW1ppW1W1PtPOI-dM1pp11PtPOI*three2npp1n1nW1ppW1W1PtPOI-dM0pp12PtPOI*three1npp1n2nW0ppW1W2PtPOI-2.*dM0pp12PtPOI*two1npp1nW1W2PtPOI-3.*dM2pp1PtPOI*two1npp1nW2ppW1PtPOI-2.*dM1pp2PtPOI-dM1pp2PtPOI*two2npp2nW1ppW2PtPOI-dM0pp3PtPOI*two1npp1nW3PtPOI-dS13mPrimePrimePtPOI+(pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePtPOI*dQnkX[0][0]+qyPrimePtPOI*dQnkY[0][0])-2.*dSnk[0][1]* (qxPrimePtPOI*dQnkX[0][0]+qyPrimePtPOI*dQnkY[0][0])+2.*(qxPrimePtPOI*dQnkX[0][2]+qyPrimePtPOI*dQnkY[0][2])-qxPrimePtPOI*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimePtPOI*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/(dM0pp111PtPOI+dM0p111PtPOI);
1417  
1418    f4WPerPtBin1n1n1n1nPOI->Fill(fPtMin+(bin-1)*dBinWidthPt,four1npp1n1n1nW1W1W1PtPOI,dM0pp111PtPOI+dM0p111PtPOI);
1419   } // end of if(dM0pp111PtPOI+dM0p111PtPOI)
1420  } // for(Int_t bin=1;bin<(fnBinsPt+1);bin++) // loop over pt-bins
1421  //...........................................................................................................
1422
1423  //...........................................................................................................  
1424  // PrimePrime Eta POI
1425  Double_t qxPrimePrimeEtaPOI=0.,qyPrimePrimeEtaPOI=0.,q2xPrimePrimeEtaPOIHere=0.,q2yPrimePrimeEtaPOIHere=0.;//add comments for these variable
1426  Double_t qxW2PrimePrimeEtaPOI=0.,qyW2PrimePrimeEtaPOI=0.,q2xW1PrimePrimeEtaPOI=0.,q2yW1PrimePrimeEtaPOI=0.;//add comments for these variable
1427  Double_t dS11mPrimePrimeEtaPOI=0.; // to be improved (name)
1428  Double_t dS12mPrimePrimeEtaPOI=0.; // to be improved (name)
1429  Double_t dS13mPrimePrimeEtaPOI=0.; // to be improved (name)
1430  Double_t mPrimePrimeEtaPOIHere=0.; // to be improved (name)
1431
1432  Double_t dM1pp11EtaPOI=0.; // to be improved (name)
1433  Double_t dM0pp111EtaPOI=0.; // to be improved (name)
1434  Double_t dM0pp12EtaPOI=0.;
1435  Double_t dM2pp1EtaPOI=0.;
1436  Double_t dM1pp2EtaPOI=0.;
1437  Double_t dM0pp3EtaPOI=0.;
1438  
1439  // Prime Eta POI
1440  Double_t qxPrimeEtaPOI=0.,qyPrimeEtaPOI=0.;
1441  Double_t mPrimeEtaPOI=0.;
1442  Double_t dM0p111EtaPOI=0.; // to be improved (name)
1443  
1444  for(Int_t bin=1;bin<(fnBinsEta+1);bin++) // loop over eta-bins 
1445  {     
1446   // q'':                    
1447   qxPrimePrimeEtaPOI = (etaReq1nPrimePrime->GetBinContent(bin))*(etaReq1nPrimePrime->GetBinEntries(bin));
1448   qyPrimePrimeEtaPOI = (etaImq1nPrimePrime->GetBinContent(bin))*(etaImq1nPrimePrime->GetBinEntries(bin)); 
1449   q2xPrimePrimeEtaPOIHere = (etaReq2nPrimePrime->GetBinContent(bin))*(etaReq2nPrimePrime->GetBinEntries(bin));  
1450   q2yPrimePrimeEtaPOIHere = (etaImq2nPrimePrime->GetBinContent(bin))*(etaImq2nPrimePrime->GetBinEntries(bin)); 
1451   
1452   qxW2PrimePrimeEtaPOI = (req1nW2PrimePrimeEta->GetBinContent(bin))*(req1nW2PrimePrimeEta->GetBinEntries(bin));
1453   qyW2PrimePrimeEtaPOI = (imq1nW2PrimePrimeEta->GetBinContent(bin))*(imq1nW2PrimePrimeEta->GetBinEntries(bin));
1454   
1455   q2xW1PrimePrimeEtaPOI = (req2nW1PrimePrimeEta->GetBinContent(bin))*(req2nW1PrimePrimeEta->GetBinEntries(bin));
1456   q2yW1PrimePrimeEtaPOI = (imq2nW1PrimePrimeEta->GetBinContent(bin))*(imq2nW1PrimePrimeEta->GetBinEntries(bin));
1457   
1458   dS11mPrimePrimeEtaPOI = (sumOfW1upTomPrimePrimeEta->GetBinContent(bin))*(sumOfW1upTomPrimePrimeEta->GetBinEntries(bin));
1459   dS12mPrimePrimeEtaPOI = (sumOfW2upTomPrimePrimeEta->GetBinContent(bin))*(sumOfW2upTomPrimePrimeEta->GetBinEntries(bin));
1460   dS13mPrimePrimeEtaPOI = (sumOfW3upTomPrimePrimeEta->GetBinContent(bin))*(sumOfW3upTomPrimePrimeEta->GetBinEntries(bin));
1461
1462   mPrimePrimeEtaPOIHere = sumOfW1upTomPrimePrimeEta->GetBinEntries(bin); // to be improved
1463       
1464   dM1pp11EtaPOI=dS11mPrimePrimeEtaPOI*(dSnk[1][0]-dSnk[0][1])-2.*dS12mPrimePrimeEtaPOI*dSnk[0][0]+2.*dS13mPrimePrimeEtaPOI;
1465   dM1pp2EtaPOI=dS11mPrimePrimeEtaPOI*dSnk[0][1]-dS13mPrimePrimeEtaPOI;
1466   dM2pp1EtaPOI=dS12mPrimePrimeEtaPOI*dSnk[0][0]-dS13mPrimePrimeEtaPOI;
1467   dM0pp3EtaPOI=mPrimePrimeEtaPOIHere*dSnk[0][2]-dS13mPrimePrimeEtaPOI;
1468   dM0pp12EtaPOI=mPrimePrimeEtaPOIHere*dSnk[0][0]*dSnk[0][1]-dM2pp1EtaPOI-dM1pp2EtaPOI-dM0pp3EtaPOI-dS13mPrimePrimeEtaPOI;
1469   dM0pp111EtaPOI=mPrimePrimeEtaPOIHere*dSnk[2][0]-3.*dM1pp11EtaPOI-3.*dM0pp12EtaPOI-3.*dM2pp1EtaPOI-3.*dM1pp2EtaPOI-dM0pp3EtaPOI-dS13mPrimePrimeEtaPOI;
1470   
1471   // q':
1472   qxPrimeEtaPOI = (etaReq1nPrime->GetBinContent(bin))*(etaReq1nPrime->GetBinEntries(bin));
1473   qyPrimeEtaPOI = (etaImq1nPrime->GetBinContent(bin))*(etaImq1nPrime->GetBinEntries(bin));
1474   
1475   mPrimeEtaPOI = etaReq1nPrime->GetBinEntries(bin); // to be improved
1476   dM0p111EtaPOI=mPrimeEtaPOI*(dSnk[2][0]-3.*dSnk[0][1]*dSnk[0][0]+2.*dSnk[0][2]);
1477  
1478   // 2-p the needed one
1479   Double_t two1n1nWPerEtaBinPOI=0.; 
1480   if((mPrimePrimeEtaPOIHere+mPrimeEtaPOI)*dSnk[0][0]-dS11mPrimePrimeEtaPOI>0)
1481   {
1482    two1n1nWPerEtaBinPOI = (qxPrimePrimeEtaPOI*dQnkX[0][0]+qyPrimePrimeEtaPOI*dQnkY[0][0]+qxPrimeEtaPOI*dQnkX[0][0]+qyPrimeEtaPOI*dQnkY[0][0]-dS11mPrimePrimeEtaPOI)/((mPrimePrimeEtaPOIHere+mPrimeEtaPOI)*dSnk[0][0]-dS11mPrimePrimeEtaPOI); 
1483   
1484    f2WPerEtaBin1n1nPOI->Fill(fEtaMin+(bin-1)*dBinWidthEta,two1n1nWPerEtaBinPOI,(mPrimePrimeEtaPOIHere+mPrimeEtaPOI)*dSnk[0][0]-dS11mPrimePrimeEtaPOI); // <2'>_{n|n} 
1485   }
1486  
1487   // 2-p the temporary one
1488   Double_t two1n1nW1ppW1W1EtaPOI=0.; // <w1 w2 w3 cos(n(phi2-phi3))> // OK!!!
1489   if(dM1pp11EtaPOI)
1490   {
1491    two1n1nW1ppW1W1EtaPOI = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*dS11mPrimePrimeEtaPOI-2.*(qxW2PrimePrimeEtaPOI*dQnkX[0][0]+qyW2PrimePrimeEtaPOI*dQnkY[0][0])-dS11mPrimePrimeEtaPOI*dSnk[0][1]+2.*dS13mPrimePrimeEtaPOI)/dM1pp11EtaPOI; // CORRECT !!! <w1 w2 w3 cos(n(phi2-phi3))>
1492   }
1493   
1494   // 2-p the temporary one
1495   Double_t two1npp1nW1W2EtaPOI=0.; // <w2 w3^2 cos(n(psi1-phi2))> // OK !!!
1496   if(dM0pp12EtaPOI)
1497   {
1498    two1npp1nW1W2EtaPOI = (dSnk[0][1]*(qxPrimePrimeEtaPOI*dQnkX[0][0]+qyPrimePrimeEtaPOI*dQnkY[0][0])-(qxW2PrimePrimeEtaPOI*dQnkX[0][0]+qyW2PrimePrimeEtaPOI*dQnkY[0][0])-dM1pp2EtaPOI-(qxPrimePrimeEtaPOI*dQnkX[0][2]+qyPrimePrimeEtaPOI*dQnkY[0][2])+dS13mPrimePrimeEtaPOI)/dM0pp12EtaPOI; // CORRECT !!! <w2 w3^2 cos(n(psi1-phi2))>
1499   }
1500   
1501   // 2-p the temporary one 
1502   Double_t two1npp1nW2ppW1EtaPOI=0.; // <w1^2 w2 cos(n(psi1-phi2))> // OK !!!
1503   if(dM2pp1EtaPOI)
1504   {
1505    two1npp1nW2ppW1EtaPOI = ((qxW2PrimePrimeEtaPOI*dQnkX[0][0]+qyW2PrimePrimeEtaPOI*dQnkY[0][0])-dS13mPrimePrimeEtaPOI)/dM2pp1EtaPOI; // CORRECT !!! <w1^2 w2 cos(n(psi1-phi2))> 
1506   }
1507   
1508   // 2-p the temporary one
1509   Double_t two2npp2nW1ppW2EtaPOI=0.; // <w1 w2^2 cos(2n(psi1-phi2))> OK !!!
1510   if(dM1pp2EtaPOI)
1511   {
1512    two2npp2nW1ppW2EtaPOI = ((q2xW1PrimePrimeEtaPOI*dQnkX[1][1]+q2yW1PrimePrimeEtaPOI*dQnkY[1][1])-dS13mPrimePrimeEtaPOI)/dM1pp2EtaPOI; // CORRECT !!! <w1 w2^2 cos(2n(psi1-phi2))> 
1513   }
1514   
1515   // 2-p the temporary one 
1516   Double_t two1npp1nW3EtaPOI=0.; // <w2^3 cos(n(psi1-phi2))> // OK !!!
1517   if(dM0pp3EtaPOI)
1518   {
1519    two1npp1nW3EtaPOI = (qxPrimePrimeEtaPOI*dQnkX[0][2]+qyPrimePrimeEtaPOI*dQnkY[0][2]-dS13mPrimePrimeEtaPOI)/dM0pp3EtaPOI; // CORRECT !!! <w2^3 cos(n(psi1-phi2))>
1520   }
1521    
1522   // 3-p the temporary one 
1523   Double_t three2npp1n1nW1ppW1W1EtaPOI=0.; // <w1 w2 w3 cos(n(2psi1-phi2-phi3))> // OK!!!
1524   if(dM1pp11EtaPOI)
1525   {
1526    three2npp1n1nW1ppW1W1EtaPOI = (q2xW1PrimePrimeEtaPOI*(dQnkX[0][0]*dQnkX[0][0]-dQnkY[0][0]*dQnkY[0][0])+2.*q2yW1PrimePrimeEtaPOI*dQnkX[0][0]*dQnkY[0][0]-2.*(qxW2PrimePrimeEtaPOI*dQnkX[0][0]+qyW2PrimePrimeEtaPOI*dQnkY[0][0])-(q2xW1PrimePrimeEtaPOI*dQnkX[1][1]+q2yW1PrimePrimeEtaPOI*dQnkY[1][1])+2.*dS13mPrimePrimeEtaPOI)/dM1pp11EtaPOI; // CORRECT !!! <w1 w2 w3 cos(n(2psi1-phi2-phi3))> 
1527   }
1528   
1529   // 3-p the temporary one 
1530   Double_t three1npp1n2nW0ppW1W2EtaPOI=0.; // <w2 w3^2 cos(n(psi1+phi2-2*phi3))> // OK!!!
1531   if(dM0pp12EtaPOI)
1532   {
1533    three1npp1n2nW0ppW1W2EtaPOI = (qxPrimePrimeEtaPOI*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])-qyPrimePrimeEtaPOI*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1])-(qxW2PrimePrimeEtaPOI*dQnkX[0][0]+qyW2PrimePrimeEtaPOI*dQnkY[0][0])-(q2xW1PrimePrimeEtaPOI*dQnkX[1][1]+q2yW1PrimePrimeEtaPOI*dQnkY[1][1])-(qxPrimePrimeEtaPOI*dQnkX[0][2]+qyPrimePrimeEtaPOI*dQnkY[0][2])+2.*dS13mPrimePrimeEtaPOI)/dM0pp12EtaPOI; // CORRECT !!! <w2 w3^2 cos(n(psi1+phi2-2.*phi3))>
1534   }
1535  
1536   /*
1537   // 4-p RP part
1538   Double_t four1npp1n1n1nW1W1W1=0.; // <w1 w2 w3 cos(n(psi1+phi1-phi2-phi3))>
1539   if(dM0pp111EtaPOI)
1540   {   
1541    four1npp1n1n1nW1W1W1 = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePrimeEtaPOI*dQnkX[0][0]+qyPrimePrimeEtaPOI*dQnkY[0][0])-2.*dM1pp11EtaPOI*two1n1nW1ppW1W1EtaPOI-dM1pp11EtaPOI*three2npp1n1nW1ppW1W1EtaPOI-dM0pp12EtaPOI*three1npp1n2nW0ppW1W2EtaPOI-2.*dM0pp12EtaPOI*two1npp1nW1W2EtaPOI-3.*dM2pp1EtaPOI*two1npp1nW2ppW1EtaPOI-2.*dM1pp2EtaPOI-dM1pp2EtaPOI*two2npp2nW1ppW2EtaPOI-dM0pp3EtaPOI*two1npp1nW3EtaPOI-dS13mPrimePrimeEtaPOI)/(dM0pp111EtaPOI); 
1542   }
1543   */
1544   /*
1545   // 4-p POI part
1546   Double_t four1npp1n1n1nW1W1W1POI=0.;
1547   if(dM0p111EtaPOI>0&&mPrimeEtaPOI>0&&nRP>0)
1548   {
1549    four1npp1n1n1nW1W1W1POI = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimeEtaPOI*dQnkX[0][0]+qyPrimeEtaPOI*dQnkY[0][0])-2.*dSnk[0][1]* (qxPrimeEtaPOI*dQnkX[0][0]+qyPrimeEtaPOI*dQnkY[0][0])+2.*(qxPrimeEtaPOI*dQnkX[0][2]+qyPrimeEtaPOI*dQnkY[0][2])-qxPrimeEtaPOI*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimeEtaPOI*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/dM0p111EtaPOI;
1550   }
1551   */
1552  
1553   // 4-p RP and POI in all combinations (full, partial and no overlap)
1554   Double_t four1npp1n1n1nW1W1W1EtaPOI=0.;
1555  
1556   if(dM0pp111EtaPOI+dM0p111EtaPOI)
1557   {
1558    four1npp1n1n1nW1W1W1EtaPOI = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePrimeEtaPOI*dQnkX[0][0]+qyPrimePrimeEtaPOI*dQnkY[0][0])-2.*dM1pp11EtaPOI*two1n1nW1ppW1W1EtaPOI-dM1pp11EtaPOI*three2npp1n1nW1ppW1W1EtaPOI-dM0pp12EtaPOI*three1npp1n2nW0ppW1W2EtaPOI-2.*dM0pp12EtaPOI*two1npp1nW1W2EtaPOI-3.*dM2pp1EtaPOI*two1npp1nW2ppW1EtaPOI-2.*dM1pp2EtaPOI-dM1pp2EtaPOI*two2npp2nW1ppW2EtaPOI-dM0pp3EtaPOI*two1npp1nW3EtaPOI-dS13mPrimePrimeEtaPOI+(pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimeEtaPOI*dQnkX[0][0]+qyPrimeEtaPOI*dQnkY[0][0])-2.*dSnk[0][1]*(qxPrimeEtaPOI*dQnkX[0][0]+qyPrimeEtaPOI*dQnkY[0][0])+2.*(qxPrimeEtaPOI*dQnkX[0][2]+qyPrimeEtaPOI*dQnkY[0][2])-qxPrimeEtaPOI*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimeEtaPOI*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/(dM0pp111EtaPOI+dM0p111EtaPOI);
1559    
1560    f4WPerEtaBin1n1n1n1nPOI->Fill(fEtaMin+(bin-1)*dBinWidthEta,four1npp1n1n1nW1W1W1EtaPOI,dM0pp111EtaPOI+dM0p111EtaPOI);
1561   }
1562  }
1563  //...........................................................................................................    
1564  
1565  
1566  
1567  
1568  
1569  
1570  
1571  
1572  
1573  
1574  
1575  
1576  
1577  
1578  // RPs
1579  ptReq1nPrime->Reset(); // to be improved 
1580  ptImq1nPrime->Reset(); // to be improved
1581  ptReq2nPrime->Reset(); // to be improved 
1582  ptImq2nPrime->Reset(); // to be improved
1583  
1584  etaReq1nPrime->Reset(); // to be improved 
1585  etaImq1nPrime->Reset(); // to be improved
1586  etaReq2nPrime->Reset(); // to be improved 
1587  etaImq2nPrime->Reset(); // to be improved
1588
1589  
1590  //...........................................................................................................
1591  // PrimePrime Pt RP
1592  Double_t qxPtRP=0.,qyPtRP=0.,q2xPtRP=0.,q2yPtRP=0.;//add comments for these variable
1593  Double_t qxW2PtRP=0.,qyW2PtRP=0.,q2xW1PtRP=0.,q2yW1PtRP=0.;//add comments for these variable
1594  Double_t dS11mPtRP=0.; // to be improved (name)
1595  Double_t dS12mPtRP=0.; // to be improved (name)
1596  Double_t dS13mPtRP=0.; // to be improved (name)
1597  Double_t mPtRP=0.; // to be improved (name)
1598
1599  Double_t dM1pp11PtRP=0.; // to be improved (name)
1600  Double_t dM0pp111PtRP=0.; // to be improved (name)
1601  Double_t dM0pp12PtRP=0.;
1602  Double_t dM2pp1PtRP=0.;
1603  Double_t dM1pp2PtRP=0.;
1604  Double_t dM0pp3PtRP=0.;
1605  
1606  // Prime Pt RP
1607  Double_t qxPrimePtRP=0.,qyPrimePtRP=0.;
1608  Double_t mPrimePtRP=0.;
1609  Double_t dM0p111PtRP=0.; // to be improved (name)
1610   
1611  for(Int_t bin=1;bin<(fnBinsPt+1);bin++) // loop over pt-bins 
1612  {     
1613   // q'':                    
1614   qxPtRP = (ptReq1n->GetBinContent(bin))*(ptReq1n->GetBinEntries(bin));
1615   qyPtRP = (ptImq1n->GetBinContent(bin))*(ptImq1n->GetBinEntries(bin)); 
1616   q2xPtRP = (ptReq2n->GetBinContent(bin))*(ptReq2n->GetBinEntries(bin));  
1617   q2yPtRP = (ptImq2n->GetBinContent(bin))*(ptImq2n->GetBinEntries(bin)); 
1618   
1619   qxW2PtRP = (req1nW2Pt->GetBinContent(bin))*(req1nW2Pt->GetBinEntries(bin));
1620   qyW2PtRP = (imq1nW2Pt->GetBinContent(bin))*(imq1nW2Pt->GetBinEntries(bin));
1621   
1622   q2xW1PtRP = (req2nW1Pt->GetBinContent(bin))*(req2nW1Pt->GetBinEntries(bin));
1623   q2yW1PtRP = (imq2nW1Pt->GetBinContent(bin))*(imq2nW1Pt->GetBinEntries(bin));
1624   
1625   dS11mPtRP = (sumOfW1upTomPt->GetBinContent(bin))*(sumOfW1upTomPt->GetBinEntries(bin));
1626   dS12mPtRP = (sumOfW2upTomPt->GetBinContent(bin))*(sumOfW2upTomPt->GetBinEntries(bin));
1627   dS13mPtRP = (sumOfW3upTomPt->GetBinContent(bin))*(sumOfW3upTomPt->GetBinEntries(bin));
1628
1629   mPtRP = sumOfW1upTomPt->GetBinEntries(bin); // to be improved
1630       
1631   dM1pp11PtRP=dS11mPtRP*(dSnk[1][0]-dSnk[0][1])-2.*dS12mPtRP*dSnk[0][0]+2.*dS13mPtRP;
1632   dM1pp2PtRP=dS11mPtRP*dSnk[0][1]-dS13mPtRP;
1633   dM2pp1PtRP=dS12mPtRP*dSnk[0][0]-dS13mPtRP;
1634   dM0pp3PtRP=mPtRP*dSnk[0][2]-dS13mPtRP;
1635   dM0pp12PtRP=mPtRP*dSnk[0][0]*dSnk[0][1]-dM2pp1PtRP-dM1pp2PtRP-dM0pp3PtRP-dS13mPtRP;
1636   dM0pp111PtRP=mPtRP*dSnk[2][0]-3.*dM1pp11PtRP-3.*dM0pp12PtRP-3.*dM2pp1PtRP-3.*dM1pp2PtRP-dM0pp3PtRP-dS13mPtRP;
1637   
1638   // q':
1639   qxPrimePtRP = (ptReq1nPrime->GetBinContent(bin))*(ptReq1nPrime->GetBinEntries(bin));
1640   qyPrimePtRP = (ptImq1nPrime->GetBinContent(bin))*(ptImq1nPrime->GetBinEntries(bin));
1641   
1642   mPrimePtRP = ptReq1nPrime->GetBinEntries(bin); // to be improved
1643   dM0p111PtRP=mPrimePtRP*(dSnk[2][0]-3.*dSnk[0][1]*dSnk[0][0]+2.*dSnk[0][2]);
1644   
1645   // 2-p the needed one
1646   Double_t two1n1nWPerPtBinRP=0.; 
1647   if((mPtRP+mPrimePtRP)*dSnk[0][0]-dS11mPtRP>0)
1648   {
1649    two1n1nWPerPtBinRP = (qxPtRP*dQnkX[0][0]+qyPtRP*dQnkY[0][0]+qxPrimePtRP*dQnkX[0][0]+qyPrimePtRP*dQnkY[0][0]-dS11mPtRP)/((mPtRP+mPrimePtRP)*dSnk[0][0]-dS11mPtRP); 
1650    f2WPerPtBin1n1nRP->Fill(fPtMin+(bin-1)*dBinWidthPt,two1n1nWPerPtBinRP,(mPtRP+mPrimePtRP)*dSnk[0][0]-dS11mPtRP);  
1651   }
1652  
1653   // 2-p temporary one
1654   Double_t two1n1nW1ppW1W1PtRP=0.; // <w1 w2 w3 cos(n(phi2-phi3))> // OK!!!
1655   if(dM1pp11PtRP)
1656   {
1657    two1n1nW1ppW1W1PtRP = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*dS11mPtRP-2.*(qxW2PtRP*dQnkX[0][0]+qyW2PtRP*dQnkY[0][0])-dS11mPtRP*dSnk[0][1]+2.*dS13mPtRP)/dM1pp11PtRP; // CORRECT !!! <w1 w2 w3 cos(n(phi2-phi3))>
1658   }
1659   
1660   // 2-p temporary one
1661   Double_t two1npp1nW1W2PtRP=0.; // <w2 w3^2 cos(n(psi1-phi2))> // OK !!!
1662   if(dM0pp12PtRP)
1663   {
1664    two1npp1nW1W2PtRP = (dSnk[0][1]*(qxPtRP*dQnkX[0][0]+qyPtRP*dQnkY[0][0])-(qxW2PtRP*dQnkX[0][0]+qyW2PtRP*dQnkY[0][0])-dM1pp2PtRP-(qxPtRP*dQnkX[0][2]+qyPtRP*dQnkY[0][2])+dS13mPtRP)/dM0pp12PtRP; // CORRECT !!! <w2 w3^2 cos(n(psi1-phi2))>
1665   }
1666   
1667   // 2-p temporary one
1668   Double_t two1npp1nW2ppW1PtRP=0.; // <w1^2 w2 cos(n(psi1-phi2))> // OK !!!
1669   if(dM2pp1PtRP)
1670   {
1671    two1npp1nW2ppW1PtRP = ((qxW2PtRP*dQnkX[0][0]+qyW2PtRP*dQnkY[0][0])-dS13mPtRP)/dM2pp1PtRP; // CORRECT !!! <w1^2 w2 cos(n(psi1-phi2))> 
1672   }
1673   
1674   // 2-p temporary one
1675   Double_t two2npp2nW1ppW2PtRP=0.; // <w1 w2^2 cos(2n(psi1-phi2))> OK !!!
1676   if(dM1pp2PtRP)
1677   {
1678    two2npp2nW1ppW2PtRP = ((q2xW1PtRP*dQnkX[1][1]+q2yW1PtRP*dQnkY[1][1])-dS13mPtRP)/dM1pp2PtRP; // CORRECT !!! <w1 w2^2 cos(2n(psi1-phi2))> 
1679   }
1680   
1681   // 2-p temporary one
1682   Double_t two1npp1nW3PtRP=0.; // <w2^3 cos(n(psi1-phi2))> // OK !!!
1683   if(dM0pp3PtRP)
1684   {
1685    two1npp1nW3PtRP = (qxPtRP*dQnkX[0][2]+qyPtRP*dQnkY[0][2]-dS13mPtRP)/dM0pp3PtRP; // CORRECT !!! <w2^3 cos(n(psi1-phi2))>
1686   }
1687    
1688   // 3-p temporary one
1689   Double_t three2npp1n1nW1ppW1W1PtRP=0.; // <w1 w2 w3 cos(n(2psi1-phi2-phi3))> // OK!!!
1690   if(dM1pp11PtRP)
1691   {
1692    three2npp1n1nW1ppW1W1PtRP = (q2xW1PtRP*(dQnkX[0][0]*dQnkX[0][0]-dQnkY[0][0]*dQnkY[0][0])+2.*q2yW1PtRP*dQnkX[0][0]*dQnkY[0][0]-2.*(qxW2PtRP*dQnkX[0][0]+qyW2PtRP*dQnkY[0][0])-(q2xW1PtRP*dQnkX[1][1]+q2yW1PtRP*dQnkY[1][1])+2.*dS13mPtRP)/dM1pp11PtRP; // CORRECT !!! <w1 w2 w3 cos(n(2psi1-phi2-phi3))> 
1693   }
1694   
1695   // 3-p temporary one
1696   Double_t three1npp1n2nW0ppW1W2PtRP=0.; // <w2 w3^2 cos(n(psi1+phi2-2*phi3))> // OK!!!
1697   if(dM0pp12PtRP)
1698   {
1699    three1npp1n2nW0ppW1W2PtRP = (qxPtRP*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])-qyPtRP*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1])-(qxW2PtRP*dQnkX[0][0]+qyW2PtRP*dQnkY[0][0])-(q2xW1PtRP*dQnkX[1][1]+q2yW1PtRP*dQnkY[1][1])-(qxPtRP*dQnkX[0][2]+qyPtRP*dQnkY[0][2])+2.*dS13mPtRP)/dM0pp12PtRP; // CORRECT !!! <w2 w3^2 cos(n(psi1+phi2-2.*phi3))>
1700   }
1701  
1702   /*
1703   // 4-p RP part
1704   Double_t four1npp1n1n1nW1W1W1=0.; // <w1 w2 w3 cos(n(psi1+phi1-phi2-phi3))>
1705   if(dM0pp111PtRP)
1706   {   
1707    four1npp1n1n1nW1W1W1 = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPtRP*dQnkX[0][0]+qyPtRP*dQnkY[0][0])-2.*dM1pp11PtRP*two1n1nW1ppW1W1PtRP-dM1pp11PtRP*three2npp1n1nW1ppW1W1PtRP-dM0pp12PtRP*three1npp1n2nW0ppW1W2PtRP-2.*dM0pp12PtRP*two1npp1nW1W2PtRP-3.*dM2pp1PtRP*two1npp1nW2ppW1PtRP-2.*dM1pp2PtRP-dM1pp2PtRP*two2npp2nW1ppW2PtRP-dM0pp3PtRP*two1npp1nW3PtRP-dS13mPtRP)/(dM0pp111PtRP); 
1708   }
1709   */
1710   
1711   /*
1712   // 4-p POI part
1713   Double_t four1npp1n1n1nW1W1W1RP=0.;
1714   if(dM0p111PtRP>0&&mPrimePtRP>0&&nRP>0)
1715   {
1716    four1npp1n1n1nW1W1W1RP = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePtRP*dQnkX[0][0]+qyPrimePtRP*dQnkY[0][0])-2.*dSnk[0][1]* (qxPrimePtRP*dQnkX[0][0]+qyPrimePtRP*dQnkY[0][0])+2.*(qxPrimePtRP*dQnkX[0][2]+qyPrimePtRP*dQnkY[0][2])-qxPrimePtRP*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimePtRP*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/dM0p111PtRP;
1717   }
1718   */
1719  
1720   // 4-p RP and POI in all combinations (full, partial and no overlap)
1721   Double_t four1npp1n1n1nW1W1W1PtRP=0.;
1722   if(dM0pp111PtRP+dM0p111PtRP)
1723   {
1724   four1npp1n1n1nW1W1W1PtRP = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPtRP*dQnkX[0][0]+qyPtRP*dQnkY[0][0])-2.*dM1pp11PtRP*two1n1nW1ppW1W1PtRP-dM1pp11PtRP*three2npp1n1nW1ppW1W1PtRP-dM0pp12PtRP*three1npp1n2nW0ppW1W2PtRP-2.*dM0pp12PtRP*two1npp1nW1W2PtRP-3.*dM2pp1PtRP*two1npp1nW2ppW1PtRP-2.*dM1pp2PtRP-dM1pp2PtRP*two2npp2nW1ppW2PtRP-dM0pp3PtRP*two1npp1nW3PtRP-dS13mPtRP+(pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimePtRP*dQnkX[0][0]+qyPrimePtRP*dQnkY[0][0])-2.*dSnk[0][1]* (qxPrimePtRP*dQnkX[0][0]+qyPrimePtRP*dQnkY[0][0])+2.*(qxPrimePtRP*dQnkX[0][2]+qyPrimePtRP*dQnkY[0][2])-qxPrimePtRP*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimePtRP*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/(dM0pp111PtRP+dM0p111PtRP);
1725  
1726    f4WPerPtBin1n1n1n1nRP->Fill(fPtMin+(bin-1)*dBinWidthPt,four1npp1n1n1nW1W1W1PtRP,dM0pp111PtRP+dM0p111PtRP);
1727   } // end of if(dM0pp111PtRP+dM0p111PtRP)
1728  } // for(Int_t bin=1;bin<(fnBinsPt+1);bin++) // loop over pt-bins
1729  
1730  delete ptReq1nPrime;
1731  delete ptImq1nPrime;
1732  delete ptReq2nPrime;
1733  delete ptImq2nPrime;
1734  delete ptReq1n;
1735  delete ptImq1n;
1736  delete ptReq2n;
1737  delete ptImq2n;
1738  
1739  delete req1nW2Pt;
1740  delete imq1nW2Pt;
1741  delete req2nW1Pt;
1742  delete imq2nW1Pt;
1743  delete sumOfW1upTomPt;
1744  delete sumOfW2upTomPt;
1745  delete sumOfW3upTomPt;
1746  //...........................................................................................................
1747
1748  //...........................................................................................................  
1749  // PrimePrime Eta RP
1750  Double_t qxEtaRP=0.,qyEtaRP=0.,q2xEtaRP=0.,q2yEtaRP=0.;//add comments for these variable
1751  Double_t qxW2EtaRP=0.,qyW2EtaRP=0.,q2xW1EtaRP=0.,q2yW1EtaRP=0.;//add comments for these variable
1752  Double_t dS11mEtaRP=0.; // to be improved (name)
1753  Double_t dS12mEtaRP=0.; // to be improved (name)
1754  Double_t dS13mEtaRP=0.; // to be improved (name)
1755  Double_t mEtaRPHere=0.; // to be improved (name)
1756
1757  Double_t dM1pp11EtaRP=0.; // to be improved (name)
1758  Double_t dM0pp111EtaRP=0.; // to be improved (name)
1759  Double_t dM0pp12EtaRP=0.;
1760  Double_t dM2pp1EtaRP=0.;
1761  Double_t dM1pp2EtaRP=0.;
1762  Double_t dM0pp3EtaRP=0.;
1763  
1764  // Prime Eta RP
1765  Double_t qxPrimeEtaRPHere=0.,qyPrimeEtaRPHere=0.;
1766  Double_t mPrimeEtaRPHere=0.;
1767  Double_t dM0p111EtaRP=0.; // to be improved (name)
1768  
1769  for(Int_t bin=1;bin<(fnBinsEta+1);bin++) // loop over eta-bins 
1770  {     
1771   // q'':                    
1772   qxEtaRP = (etaReq1n->GetBinContent(bin))*(etaReq1n->GetBinEntries(bin));
1773   qyEtaRP = (etaImq1n->GetBinContent(bin))*(etaImq1n->GetBinEntries(bin)); 
1774   q2xEtaRP = (etaReq2n->GetBinContent(bin))*(etaReq2n->GetBinEntries(bin));  
1775   q2yEtaRP = (etaImq2n->GetBinContent(bin))*(etaImq2n->GetBinEntries(bin)); 
1776   
1777   qxW2EtaRP = (req1nW2Eta->GetBinContent(bin))*(req1nW2Eta->GetBinEntries(bin));
1778   qyW2EtaRP = (imq1nW2Eta->GetBinContent(bin))*(imq1nW2Eta->GetBinEntries(bin));
1779   
1780   q2xW1EtaRP = (req2nW1Eta->GetBinContent(bin))*(req2nW1Eta->GetBinEntries(bin));
1781   q2yW1EtaRP = (imq2nW1Eta->GetBinContent(bin))*(imq2nW1Eta->GetBinEntries(bin));
1782   
1783   dS11mEtaRP = (sumOfW1upTomEta->GetBinContent(bin))*(sumOfW1upTomEta->GetBinEntries(bin));
1784   dS12mEtaRP = (sumOfW2upTomEta->GetBinContent(bin))*(sumOfW2upTomEta->GetBinEntries(bin));
1785   dS13mEtaRP = (sumOfW3upTomEta->GetBinContent(bin))*(sumOfW3upTomEta->GetBinEntries(bin));
1786
1787   mEtaRPHere = sumOfW1upTomEta->GetBinEntries(bin); // to be improved
1788       
1789   dM1pp11EtaRP=dS11mEtaRP*(dSnk[1][0]-dSnk[0][1])-2.*dS12mEtaRP*dSnk[0][0]+2.*dS13mEtaRP;
1790   dM1pp2EtaRP=dS11mEtaRP*dSnk[0][1]-dS13mEtaRP;
1791   dM2pp1EtaRP=dS12mEtaRP*dSnk[0][0]-dS13mEtaRP;
1792   dM0pp3EtaRP=mEtaRPHere*dSnk[0][2]-dS13mEtaRP;
1793   dM0pp12EtaRP=mEtaRPHere*dSnk[0][0]*dSnk[0][1]-dM2pp1EtaRP-dM1pp2EtaRP-dM0pp3EtaRP-dS13mEtaRP;
1794   dM0pp111EtaRP=mEtaRPHere*dSnk[2][0]-3.*dM1pp11EtaRP-3.*dM0pp12EtaRP-3.*dM2pp1EtaRP-3.*dM1pp2EtaRP-dM0pp3EtaRP-dS13mEtaRP;
1795   
1796   // q':
1797   qxPrimeEtaRPHere = (etaReq1nPrime->GetBinContent(bin))*(etaReq1nPrime->GetBinEntries(bin));
1798   qyPrimeEtaRPHere = (etaImq1nPrime->GetBinContent(bin))*(etaImq1nPrime->GetBinEntries(bin));
1799   
1800   mPrimeEtaRPHere = etaReq1nPrime->GetBinEntries(bin); // to be improved
1801   dM0p111EtaRP=mPrimeEtaRPHere*(dSnk[2][0]-3.*dSnk[0][1]*dSnk[0][0]+2.*dSnk[0][2]);
1802  
1803   // 2-p the needed one
1804   Double_t two1n1nWPerEtaBinRP=0.; 
1805   if((mEtaRPHere+mPrimeEtaRPHere)*dSnk[0][0]-dS11mEtaRP>0)
1806   {
1807    two1n1nWPerEtaBinRP = (qxEtaRP*dQnkX[0][0]+qyEtaRP*dQnkY[0][0]+qxPrimeEtaRPHere*dQnkX[0][0]+qyPrimeEtaRPHere*dQnkY[0][0]-dS11mEtaRP)/((mEtaRPHere+mPrimeEtaRPHere)*dSnk[0][0]-dS11mEtaRP); 
1808   
1809    f2WPerEtaBin1n1nRP->Fill(fEtaMin+(bin-1)*dBinWidthEta,two1n1nWPerEtaBinRP,(mEtaRPHere+mPrimeEtaRPHere)*dSnk[0][0]-dS11mEtaRP); // <2'>_{n|n} 
1810   }
1811  
1812   // 2-p the temporary one
1813   Double_t two1n1nW1ppW1W1EtaRP=0.; // <w1 w2 w3 cos(n(phi2-phi3))> // OK!!!
1814   if(dM1pp11EtaRP)
1815   {
1816    two1n1nW1ppW1W1EtaRP = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*dS11mEtaRP-2.*(qxW2EtaRP*dQnkX[0][0]+qyW2EtaRP*dQnkY[0][0])-dS11mEtaRP*dSnk[0][1]+2.*dS13mEtaRP)/dM1pp11EtaRP; // CORRECT !!! <w1 w2 w3 cos(n(phi2-phi3))>
1817   }
1818   
1819   // 2-p the temporary one
1820   Double_t two1npp1nW1W2EtaRP=0.; // <w2 w3^2 cos(n(psi1-phi2))> // OK !!!
1821   if(dM0pp12EtaRP)
1822   {
1823    two1npp1nW1W2EtaRP = (dSnk[0][1]*(qxEtaRP*dQnkX[0][0]+qyEtaRP*dQnkY[0][0])-(qxW2EtaRP*dQnkX[0][0]+qyW2EtaRP*dQnkY[0][0])-dM1pp2EtaRP-(qxEtaRP*dQnkX[0][2]+qyEtaRP*dQnkY[0][2])+dS13mEtaRP)/dM0pp12EtaRP; // CORRECT !!! <w2 w3^2 cos(n(psi1-phi2))>
1824   }
1825   
1826   // 2-p the temporary one 
1827   Double_t two1npp1nW2ppW1EtaRP=0.; // <w1^2 w2 cos(n(psi1-phi2))> // OK !!!
1828   if(dM2pp1EtaRP)
1829   {
1830    two1npp1nW2ppW1EtaRP = ((qxW2EtaRP*dQnkX[0][0]+qyW2EtaRP*dQnkY[0][0])-dS13mEtaRP)/dM2pp1EtaRP; // CORRECT !!! <w1^2 w2 cos(n(psi1-phi2))> 
1831   }
1832   
1833   // 2-p the temporary one
1834   Double_t two2npp2nW1ppW2EtaRP=0.; // <w1 w2^2 cos(2n(psi1-phi2))> OK !!!
1835   if(dM1pp2EtaRP)
1836   {
1837    two2npp2nW1ppW2EtaRP = ((q2xW1EtaRP*dQnkX[1][1]+q2yW1EtaRP*dQnkY[1][1])-dS13mEtaRP)/dM1pp2EtaRP; // CORRECT !!! <w1 w2^2 cos(2n(psi1-phi2))> 
1838   }
1839   
1840   // 2-p the temporary one 
1841   Double_t two1npp1nW3EtaRP=0.; // <w2^3 cos(n(psi1-phi2))> // OK !!!
1842   if(dM0pp3EtaRP)
1843   {
1844    two1npp1nW3EtaRP = (qxEtaRP*dQnkX[0][2]+qyEtaRP*dQnkY[0][2]-dS13mEtaRP)/dM0pp3EtaRP; // CORRECT !!! <w2^3 cos(n(psi1-phi2))>
1845   }
1846    
1847   // 3-p the temporary one 
1848   Double_t three2npp1n1nW1ppW1W1EtaRP=0.; // <w1 w2 w3 cos(n(2psi1-phi2-phi3))> // OK!!!
1849   if(dM1pp11EtaRP)
1850   {
1851    three2npp1n1nW1ppW1W1EtaRP = (q2xW1EtaRP*(dQnkX[0][0]*dQnkX[0][0]-dQnkY[0][0]*dQnkY[0][0])+2.*q2yW1EtaRP*dQnkX[0][0]*dQnkY[0][0]-2.*(qxW2EtaRP*dQnkX[0][0]+qyW2EtaRP*dQnkY[0][0])-(q2xW1EtaRP*dQnkX[1][1]+q2yW1EtaRP*dQnkY[1][1])+2.*dS13mEtaRP)/dM1pp11EtaRP; // CORRECT !!! <w1 w2 w3 cos(n(2psi1-phi2-phi3))> 
1852   }
1853   
1854   // 3-p the temporary one 
1855   Double_t three1npp1n2nW0ppW1W2EtaRP=0.; // <w2 w3^2 cos(n(psi1+phi2-2*phi3))> // OK!!!
1856   if(dM0pp12EtaRP)
1857   {
1858    three1npp1n2nW0ppW1W2EtaRP = (qxEtaRP*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])-qyEtaRP*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1])-(qxW2EtaRP*dQnkX[0][0]+qyW2EtaRP*dQnkY[0][0])-(q2xW1EtaRP*dQnkX[1][1]+q2yW1EtaRP*dQnkY[1][1])-(qxEtaRP*dQnkX[0][2]+qyEtaRP*dQnkY[0][2])+2.*dS13mEtaRP)/dM0pp12EtaRP; // CORRECT !!! <w2 w3^2 cos(n(psi1+phi2-2.*phi3))>
1859   }
1860  
1861   /*
1862   // 4-p RP part
1863   Double_t four1npp1n1n1nW1W1W1=0.; // <w1 w2 w3 cos(n(psi1+phi1-phi2-phi3))>
1864   if(dM0pp111EtaRP)
1865   {   
1866    four1npp1n1n1nW1W1W1 = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxEtaRP*dQnkX[0][0]+qyEtaRP*dQnkY[0][0])-2.*dM1pp11EtaRP*two1n1nW1ppW1W1EtaRP-dM1pp11EtaRP*three2npp1n1nW1ppW1W1EtaRP-dM0pp12EtaRP*three1npp1n2nW0ppW1W2EtaRP-2.*dM0pp12EtaRP*two1npp1nW1W2EtaRP-3.*dM2pp1EtaRP*two1npp1nW2ppW1EtaRP-2.*dM1pp2EtaRP-dM1pp2EtaRP*two2npp2nW1ppW2EtaRP-dM0pp3EtaRP*two1npp1nW3EtaRP-dS13mEtaRP)/(dM0pp111EtaRP); 
1867   }
1868   */
1869   /*
1870   // 4-p POI part
1871   Double_t four1npp1n1n1nW1W1W1RP=0.;
1872   if(dM0p111EtaRP>0&&mPrimeEtaRPHere>0&&nRP>0)
1873   {
1874    four1npp1n1n1nW1W1W1RP = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimeEtaRPHere*dQnkX[0][0]+qyPrimeEtaRPHere*dQnkY[0][0])-2.*dSnk[0][1]* (qxPrimeEtaRPHere*dQnkX[0][0]+qyPrimeEtaRPHere*dQnkY[0][0])+2.*(qxPrimeEtaRPHere*dQnkX[0][2]+qyPrimeEtaRPHere*dQnkY[0][2])-qxPrimeEtaRPHere*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimeEtaRPHere*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/dM0p111EtaRP;
1875   }
1876   */
1877  
1878   // 4-p RP and POI in all combinations (full, partial and no overlap)
1879   Double_t four1npp1n1n1nW1W1W1EtaRP=0.;
1880  
1881   if(dM0pp111EtaRP+dM0p111EtaRP)
1882   {
1883    four1npp1n1n1nW1W1W1EtaRP = ((pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxEtaRP*dQnkX[0][0]+qyEtaRP*dQnkY[0][0])-2.*dM1pp11EtaRP*two1n1nW1ppW1W1EtaRP-dM1pp11EtaRP*three2npp1n1nW1ppW1W1EtaRP-dM0pp12EtaRP*three1npp1n2nW0ppW1W2EtaRP-2.*dM0pp12EtaRP*two1npp1nW1W2EtaRP-3.*dM2pp1EtaRP*two1npp1nW2ppW1EtaRP-2.*dM1pp2EtaRP-dM1pp2EtaRP*two2npp2nW1ppW2EtaRP-dM0pp3EtaRP*two1npp1nW3EtaRP-dS13mEtaRP+(pow(dQnkX[0][0],2.)+pow(dQnkY[0][0],2.))*(qxPrimeEtaRPHere*dQnkX[0][0]+qyPrimeEtaRPHere*dQnkY[0][0])-2.*dSnk[0][1]*(qxPrimeEtaRPHere*dQnkX[0][0]+qyPrimeEtaRPHere*dQnkY[0][0])+2.*(qxPrimeEtaRPHere*dQnkX[0][2]+qyPrimeEtaRPHere*dQnkY[0][2])-qxPrimeEtaRPHere*(dQnkX[0][0]*dQnkX[1][1]+dQnkY[0][0]*dQnkY[1][1])+qyPrimeEtaRPHere*(dQnkY[0][0]*dQnkX[1][1]-dQnkX[0][0]*dQnkY[1][1]))/(dM0pp111EtaRP+dM0p111EtaRP);
1884    
1885    f4WPerEtaBin1n1n1n1nRP->Fill(fEtaMin+(bin-1)*dBinWidthEta,four1npp1n1n1nW1W1W1EtaRP,dM0pp111EtaRP+dM0p111EtaRP);
1886   }
1887  }
1888
1889  delete etaReq1nPrime;
1890  delete etaImq1nPrime;
1891  delete etaReq2nPrime;
1892  delete etaImq2nPrime;
1893  delete etaReq1n;
1894  delete etaImq1n;
1895  delete etaReq2n;
1896  delete etaImq2n;
1897
1898  delete req1nW2Eta;
1899  delete imq1nW2Eta;
1900  delete req2nW1Eta;
1901  delete imq2nW1Eta;
1902  delete sumOfW1upTomEta;
1903  delete sumOfW2upTomEta;
1904  delete sumOfW3upTomEta;
1905  //........................................................................................................... 
1906  
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924  
1925  
1926  
1927  
1928  
1929  
1930  
1931  
1932  
1933  
1934  
1935  
1936  
1937  
1938  
1939  
1940  
1941  
1942  
1943  
1944  
1945  
1946  
1947  
1948  
1949  
1950  
1951  
1952  
1953  
1954  
1955  
1956  
1957  
1958  
1959  
1960  
1961  
1962  
1963  
1964  
1965  
1966  
1967  
1968  
1969  
1970  
1971  
1972  
1973  
1974  
1975  
1976  
1977  
1978  
1979  
1980  
1981  
1982  
1983  
1984  
1985  
1986  
1987  
1988  
1989  
1990  
1991  
1992  
1993  
1994  
1995  
1996  
1997  
1998  
1999  
2000  
2001  
2002  
2003  
2004  
2005  
2006  
2007  
2008  
2009  
2010  
2011  
2012  
2013  
2014  
2015  
2016  
2017  
2018  
2019  
2020  
2021  
2022  
2023  
2024  
2025  
2026  
2027  
2028  
2029  
2030  
2031  
2032  
2033  
2034  
2035  
2036  
2037  
2038  
2039  
2040  
2041  
2042  
2043  
2044  
2045  
2046  
2047  
2048  
2049  
2050  
2051  
2052  
2053  
2054  
2055  
2056  
2057  
2058  
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084  /*
2085
2086  //if(bNestedLoops)to be improved
2087  //{ to be improved               
2088  //-------------------------------------------------------------------------------------------------------------------------------- 
2089  //
2090  //                                          **********************
2091  //                                          **** NESTED LOOPS ****
2092  //                                          ********************** 
2093  //
2094  // Remark 1: multi-particle correlations calculated with nested loops are stored in fDirectCorrelations.
2095  // Remark 2: binning of fDirectCorrelations: bins 0..100 - correlations needed for integrated flow; bins 100..200 - correlations needed for differential flow (taking as an example bin 0.5 < pt < 0.6)
2096  //
2097  // binning details of fDirectCorrelations (integrated flow):
2098  //..........................................................................
2099  // 1st bin: weighted <2>_{n|n}   = <w1   w2   cos( n*(phi1-phi2))>
2100  // 2nd bin: weighted <2>_{2n|2n} = <w1^2 w2^2 cos(2n*(phi1-phi2))>
2101  // 3rd bin: weighted <2>_{3n|3n} = <w1^3 w2^3 cos(3n*(phi1-phi2))>
2102  // 4th bin: weighted <2>_{4n|4n} = <w1^4 w2^4 cos(4n*(phi1-phi2))>
2103  // 5th bin: weighted <2>_{n|n} = <w1^3 w2 cos(n*(phi1-phi2))>
2104  
2105  // 11th bin: weighted <3>_{2n|n,n} = <w1^2 w2 w3 cos(n*(2phi1-phi2-phi3))>
2106  //..........................................................................
2107  
2108
2109  Double_t phi1=0., phi2=0., phi3=0., phi4=0., phi5=0., phi6=0., phi7=0., phi8=0.;
2110  Double_t wPhi1=1., wPhi2=1., wPhi3=1., wPhi4=1., wPhi5=1., wPhi6=1., wPhi7=1., wPhi8=1.;
2111
2112
2113
2114  
2115
2116
2117  Double_t tempLoop = 0.;
2118  Int_t tempCounter = 0;
2119  
2120  
2121  
2122  
2123  
2124  for(Int_t i1=0;i1<dMult;i1++)
2125  {
2126   fTrack=anEvent->GetTrack(i1);
2127   phi1=fTrack->Phi();
2128   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2129   for(Int_t i2=0;i2<dMult;i2++)
2130   {
2131    if(i2==i1)continue;
2132    fTrack=anEvent->GetTrack(i2);
2133    phi2=fTrack->Phi();
2134    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2135    // 2-p
2136    fDirectCorrelations->Fill(0.,cos(n*(phi1-phi2)),wPhi1*wPhi2);                  // <w1   w2   cos( n*(phi1-phi2))>
2137    fDirectCorrelations->Fill(1.,cos(2.*n*(phi1-phi2)),pow(wPhi1,2)*pow(wPhi2,2)); // <w1^2 w2^2 cos(2n*(phi1-phi2))>
2138    fDirectCorrelations->Fill(2.,cos(3.*n*(phi1-phi2)),pow(wPhi1,3)*pow(wPhi2,3)); // <w1^3 w2^3 cos(3n*(phi1-phi2))>
2139    fDirectCorrelations->Fill(3.,cos(4.*n*(phi1-phi2)),pow(wPhi1,4)*pow(wPhi2,4)); // <w1^4 w2^4 cos(4n*(phi1-phi2))> 
2140    fDirectCorrelations->Fill(4.,cos(n*(phi1-phi2)),pow(wPhi1,3)*wPhi2); // <w1^3 w2 cos(n*(phi1-phi2))>
2141   }
2142  }  
2143  
2144
2145
2146
2147
2148
2149  for(Int_t i1=0;i1<dMult;i1++)
2150  {
2151   fTrack=anEvent->GetTrack(i1);
2152   phi1=fTrack->Phi();
2153   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2154   for(Int_t i2=0;i2<dMult;i2++)
2155   {
2156    if(i2==i1)continue;
2157    fTrack=anEvent->GetTrack(i2);
2158    phi2=fTrack->Phi();
2159    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2160    for(Int_t i3=0;i3<dMult;i3++)
2161    {
2162     if(i3==i1||i3==i2)continue;
2163     fTrack=anEvent->GetTrack(i3);
2164     phi3=fTrack->Phi();
2165     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2166     // 2-p
2167     fDirectCorrelations->Fill(5.,cos(n*(phi1-phi2)),wPhi1*wPhi2*pow(wPhi3,2)); // <w1 w2 w3^2 cos(n*(phi1-phi2))>
2168     
2169     // 3-p
2170     fDirectCorrelations->Fill(10.,cos(2.*n*phi1-n*(phi2+phi3)),pow(wPhi1,2)*wPhi2*wPhi3); // <w1^2 w2 w3 cos(n*(2phi1-phi2-phi3))>
2171     
2172     
2173     fDirectCorrelations->Fill(6.,cos(3.*n*phi1-2.*n*phi2-n*phi3),pow(wPhi1,3)*pow(wPhi2,2)*wPhi3);           //<3>_{3n|2n,n}
2174     fDirectCorrelations->Fill(7.,cos(4.*n*phi1-2.*n*phi2-2.*n*phi3),pow(wPhi1,4)*pow(wPhi2,2)*pow(wPhi3,2)); //<3>_{4n|2n,2n}
2175     fDirectCorrelations->Fill(8.,cos(4.*n*phi1-3.*n*phi2-n*phi3),pow(wPhi1,4)*pow(wPhi2,3)*wPhi3);           //<3>_{4n|3n,n}
2176     
2177     
2178    }
2179   }
2180  }
2181  
2182
2183  
2184  
2185
2186
2187   
2188  //<4>_{n,n|n,n}, <4>_{2n,n|2n,n}, <4>_{2n,2n|2n,2n}, <4>_{3n|n,n,n}, <4>_{3n,n|3n,n}, <4>_{3n,n|2n,2n} and <4>_{4n|2n,n,n} 
2189  for(Int_t i1=0;i1<dMult;i1++)
2190  {
2191   fTrack=anEvent->GetTrack(i1);
2192   phi1=fTrack->Phi();
2193   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2194   for(Int_t i2=0;i2<dMult;i2++)
2195   {
2196    if(i2==i1)continue;
2197    fTrack=anEvent->GetTrack(i2);
2198    phi2=fTrack->Phi();
2199    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2200    for(Int_t i3=0;i3<dMult;i3++)
2201    {
2202     if(i3==i1||i3==i2)continue;
2203     fTrack=anEvent->GetTrack(i3);
2204     phi3=fTrack->Phi();
2205     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2206     for(Int_t i4=0;i4<dMult;i4++)
2207     {
2208      if(i4==i1||i4==i2||i4==i3)continue;
2209      fTrack=anEvent->GetTrack(i4);
2210      phi4=fTrack->Phi();
2211      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2212      fDirectCorrelations->Fill(20.,cos(n*phi1+n*phi2-n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4); // <4>_{n,n|n,n} = <w1 w2 w3 w4 cos(n*(phi1+phi2-phi3-phi4))>
2213      
2214     
2215      fDirectCorrelations->Fill(11.,cos(2.*n*phi1+n*phi2-2.*n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);      //<4>_{2n,n|2n,n}
2216      fDirectCorrelations->Fill(12.,cos(2.*n*phi1+2*n*phi2-2.*n*phi3-2.*n*phi4),wPhi1*wPhi2*wPhi3*wPhi4); //<4>_{2n,2n|2n,2n}
2217      fDirectCorrelations->Fill(13.,cos(3.*n*phi1-n*phi2-n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);         //<4>_{3n|n,n,n}
2218      fDirectCorrelations->Fill(14.,cos(3.*n*phi1+n*phi2-3.*n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);      //<4>_{3n,n|3n,n}   
2219      fDirectCorrelations->Fill(15.,cos(3.*n*phi1+n*phi2-2.*n*phi3-2.*n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);   //<4>_{3n,n|2n,2n}
2220      fDirectCorrelations->Fill(16.,cos(4.*n*phi1-2.*n*phi2-n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);      //<4>_{4n|2n,n,n}
2221     
2222     }  
2223    }
2224   }
2225  }
2226  
2227  
2228  
2229
2230  
2231  //<5>_{2n,n,n,n,n}, //<5>_{2n,2n|2n,n,n}, <5>_{3n,n|2n,n,n} and <5>_{4n|n,n,n,n}
2232  for(Int_t i1=0;i1<dMult;i1++)
2233  {
2234   //cout<<"i1 = "<<i1<<endl;
2235   fTrack=anEvent->GetTrack(i1);
2236   phi1=fTrack->Phi();
2237   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2238   for(Int_t i2=0;i2<dMult;i2++)
2239   {
2240    if(i2==i1)continue;
2241    fTrack=anEvent->GetTrack(i2);
2242    phi2=fTrack->Phi();
2243    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2244    for(Int_t i3=0;i3<dMult;i3++)
2245    {
2246     if(i3==i1||i3==i2)continue;
2247     fTrack=anEvent->GetTrack(i3);
2248     phi3=fTrack->Phi();
2249     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2250     for(Int_t i4=0;i4<dMult;i4++)
2251     {
2252      if(i4==i1||i4==i2||i4==i3)continue;
2253      fTrack=anEvent->GetTrack(i4);
2254      phi4=fTrack->Phi();
2255      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2256      for(Int_t i5=0;i5<dMult;i5++)
2257      {
2258       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2259       fTrack=anEvent->GetTrack(i5);
2260       phi5=fTrack->Phi();
2261       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2262       fDirectCorrelations->Fill(18.,cos(2.*n*phi1+n*phi2-n*phi3-n*phi4-n*phi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5);       //<5>_{2n,n|n,n,n}
2263       fDirectCorrelations->Fill(19.,cos(2.*n*phi1+2.*n*phi2-2.*n*phi3-n*phi4-n*phi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5); //<5>_{2n,2n|2n,n,n}
2264       fDirectCorrelations->Fill(20.,cos(3.*n*phi1+n*phi2-2.*n*phi3-n*phi4-n*phi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5);    //<5>_{3n,n|2n,n,n}
2265       fDirectCorrelations->Fill(21.,cos(4.*n*phi1-n*phi2-n*phi3-n*phi4-n*phi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5);       //<5>_{4n|n,n,n,n}
2266      }
2267     }  
2268    }
2269   }
2270  }
2271  
2272  
2273  //<6>_{n,n,n,n,n,n}, <6>_{2n,n,n|2n,n,n}, <6>_{2n,2n|n,n,n,n} and <6>_{3n,n|n,n,n,n}
2274  for(Int_t i1=0;i1<dMult;i1++)
2275  {
2276   //cout<<"i1 = "<<i1<<endl;
2277   fTrack=anEvent->GetTrack(i1);
2278   phi1=fTrack->Phi();
2279   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2280   for(Int_t i2=0;i2<dMult;i2++)
2281   {
2282    if(i2==i1)continue;
2283    fTrack=anEvent->GetTrack(i2);
2284    phi2=fTrack->Phi();
2285    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2286    for(Int_t i3=0;i3<dMult;i3++)
2287    {
2288     if(i3==i1||i3==i2)continue;
2289     fTrack=anEvent->GetTrack(i3);
2290     phi3=fTrack->Phi();
2291     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2292     for(Int_t i4=0;i4<dMult;i4++)
2293     {
2294      if(i4==i1||i4==i2||i4==i3)continue;
2295      fTrack=anEvent->GetTrack(i4);
2296      phi4=fTrack->Phi();
2297      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2298      for(Int_t i5=0;i5<dMult;i5++)
2299      {
2300       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2301       fTrack=anEvent->GetTrack(i5);
2302       phi5=fTrack->Phi();
2303       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2304       for(Int_t i6=0;i6<dMult;i6++)
2305       {
2306        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2307        fTrack=anEvent->GetTrack(i6);
2308        phi6=fTrack->Phi(); 
2309        if(phiWeights) wPhi6 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi6*nBinsPhi/TMath::TwoPi())));
2310        fDirectCorrelations->Fill(23.,cos(n*phi1+n*phi2+n*phi3-n*phi4-n*phi5-n*phi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6);       //<6>_{n,n,n|n,n,n}
2311        fDirectCorrelations->Fill(24.,cos(2.*n*phi1+n*phi2+n*phi3-2.*n*phi4-n*phi5-n*phi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6); //<6>_{2n,n,n|2n,n,n}
2312        fDirectCorrelations->Fill(25.,cos(2.*n*phi1+2.*n*phi2-n*phi3-n*phi4-n*phi5-n*phi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6); //<6>_{2n,2n|n,n,n,n}
2313        fDirectCorrelations->Fill(26.,cos(3.*n*phi1+n*phi2-n*phi3-n*phi4-n*phi5-n*phi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6);    //<6>_{3n,n|n,n,n,n}     
2314       } 
2315      }
2316     }  
2317    }
2318   }
2319  }
2320
2321
2322  
2323  //<7>_{2n,n,n|n,n,n,n}
2324  for(Int_t i1=0;i1<dMult;i1++)
2325  {
2326   //cout<<"i1 = "<<i1<<endl;
2327   fTrack=anEvent->GetTrack(i1);
2328   phi1=fTrack->Phi();
2329   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2330   for(Int_t i2=0;i2<dMult;i2++)
2331   {
2332    if(i2==i1)continue;
2333    fTrack=anEvent->GetTrack(i2);
2334    phi2=fTrack->Phi();
2335    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2336    for(Int_t i3=0;i3<dMult;i3++)
2337    {
2338     if(i3==i1||i3==i2)continue;
2339     fTrack=anEvent->GetTrack(i3);
2340     phi3=fTrack->Phi();
2341     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2342     for(Int_t i4=0;i4<dMult;i4++)
2343     {
2344      if(i4==i1||i4==i2||i4==i3)continue;
2345      fTrack=anEvent->GetTrack(i4);
2346      phi4=fTrack->Phi();
2347      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2348      for(Int_t i5=0;i5<dMult;i5++)
2349      {
2350       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2351       fTrack=anEvent->GetTrack(i5);
2352       phi5=fTrack->Phi();
2353       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2354       for(Int_t i6=0;i6<dMult;i6++)
2355       {
2356        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2357        fTrack=anEvent->GetTrack(i6);
2358        phi6=fTrack->Phi(); 
2359        if(phiWeights) wPhi6 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi6*nBinsPhi/TMath::TwoPi())));
2360        for(Int_t i7=0;i7<dMult;i7++)
2361        {
2362         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2363         fTrack=anEvent->GetTrack(i7);
2364         phi7=fTrack->Phi(); 
2365         if(phiWeights) wPhi7 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi7*nBinsPhi/TMath::TwoPi())));
2366         fDirectCorrelations->Fill(28.,cos(2.*n*phi1+n*phi2+n*phi3-n*phi4-n*phi5-n*phi6-n*phi7),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6*wPhi7);//<7>_{2n,n,n|n,n,n,n}
2367        } 
2368       } 
2369      }
2370     }  
2371    }
2372   }
2373  }
2374  
2375  
2376  
2377  //<8>_{n,n,n,n|n,n,n,n}
2378  for(Int_t i1=0;i1<dMult;i1++)
2379  {
2380   cout<<"i1 = "<<i1<<endl;
2381   fTrack=anEvent->GetTrack(i1);
2382   phi1=fTrack->Phi();
2383   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2384   for(Int_t i2=0;i2<dMult;i2++)
2385   {
2386    if(i2==i1)continue;
2387    fTrack=anEvent->GetTrack(i2);
2388    phi2=fTrack->Phi();
2389    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2390    for(Int_t i3=0;i3<dMult;i3++)
2391    {
2392     if(i3==i1||i3==i2)continue;
2393     fTrack=anEvent->GetTrack(i3);
2394     phi3=fTrack->Phi();
2395     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2396     for(Int_t i4=0;i4<dMult;i4++)
2397     {
2398      if(i4==i1||i4==i2||i4==i3)continue;
2399      fTrack=anEvent->GetTrack(i4);
2400      phi4=fTrack->Phi();
2401      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2402      for(Int_t i5=0;i5<dMult;i5++)
2403      {
2404       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2405       fTrack=anEvent->GetTrack(i5);
2406       phi5=fTrack->Phi();
2407       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2408       for(Int_t i6=0;i6<dMult;i6++)
2409       {
2410        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2411        fTrack=anEvent->GetTrack(i6);
2412        phi6=fTrack->Phi();
2413        if(phiWeights) wPhi6 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi6*nBinsPhi/TMath::TwoPi()))); 
2414        for(Int_t i7=0;i7<dMult;i7++)
2415        {
2416         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2417         fTrack=anEvent->GetTrack(i7);
2418         phi7=fTrack->Phi();
2419         if(phiWeights) wPhi7 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi7*nBinsPhi/TMath::TwoPi()))); 
2420         for(Int_t i8=0;i8<dMult;i8++)
2421         {
2422          if(i8==i1||i8==i2||i8==i3||i8==i4||i8==i5||i8==i6||i8==i7)continue;
2423          fTrack=anEvent->GetTrack(i8);
2424          phi8=fTrack->Phi();
2425          if(phiWeights) wPhi8 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi8*nBinsPhi/TMath::TwoPi())));  
2426          fDirectCorrelations->Fill(30.,cos(n*phi1+n*phi2+n*phi3+n*phi4-n*phi5-n*phi6-n*phi7-n*phi8),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6*wPhi7*wPhi8);//<8>_{n,n,n,n|n,n,n,n}
2427         } 
2428        } 
2429       } 
2430      }
2431     }  
2432    }
2433   }
2434  }
2435  
2436  */
2437
2438  
2439
2440  // binning details of fDirectCorrelations (differential flow):
2441  //
2442  //101st bin: <2'>_{n|n}
2443  
2444  
2445  /*
2446  
2447  //<2'>_{n|n}
2448  for(Int_t i1=0;i1<nPrim;i1++)
2449  {
2450   fTrack=anEvent->GetTrack(i1);
2451   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2452   phi1=fTrack->Phi();
2453   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2454   for(Int_t i2=0;i2<nPrim;i2++)
2455   {
2456    if(i2==i1)continue;
2457    fTrack=anEvent->GetTrack(i2);
2458    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition
2459    phi2=fTrack->Phi();
2460    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2461     //cout<<"1st = "<<i1<<"     "<< (anEvent->GetTrack(i1))->Eta() << " " << (anEvent->GetTrack(i1))->Pt()<<endl;
2462     //cout<<"2nd = "<<i2<<"     "<< (anEvent->GetTrack(i2))->Eta() << " " << (anEvent->GetTrack(i2))->Pt()<<endl; 
2463    //fill the fDirectCorrelations:    
2464    fDirectCorrelations->Fill(100.,cos(1.*n*(phi1-phi2)),wPhi2);//<2'>_{n,n}
2465    fDirectCorrelations->Fill(103.,cos(1.*n*(phi1-phi2)),pow(wPhi1,2)*wPhi2);//<2'>_{n,n}
2466    fDirectCorrelations->Fill(104.,cos(2.*n*(phi1-phi2)),wPhi1*pow(wPhi2,2));//<2'>_{n,n}
2467    fDirectCorrelations->Fill(105.,cos(1.*n*(phi1-phi2)),pow(wPhi2,3));//<2'>_{n,n}
2468    
2469    
2470    
2471    fDirectCorrelations->Fill(41.,cos(2.*n*(phi1-phi2)),1);//<2'>_{2n,2n}
2472    fDirectCorrelations->Fill(42.,cos(3.*n*(phi1-phi2)),1);//<2'>_{3n,3n}
2473    fDirectCorrelations->Fill(43.,cos(4.*n*(phi1-phi2)),1);//<2'>_{4n,4n}   
2474     
2475   }//end of for(Int_t i2=0;i2<nPrim;i2++)
2476  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2477
2478  */ 
2479   
2480  
2481  
2482  /*
2483  
2484  //<3'>_{2n|n,n}
2485  for(Int_t i1=0;i1<nPrim;i1++)
2486  {
2487   fTrack=anEvent->GetTrack(i1);
2488   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2489   phi1=fTrack->Phi();
2490   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2491   for(Int_t i2=0;i2<nPrim;i2++)
2492   {
2493    if(i2==i1)continue;
2494    fTrack=anEvent->GetTrack(i2);
2495    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition
2496    phi2=fTrack->Phi();
2497    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2498    for(Int_t i3=0;i3<nPrim;i3++)
2499    {
2500     if(i3==i1||i3==i2)continue;
2501     fTrack=anEvent->GetTrack(i3);
2502     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition
2503     phi3=fTrack->Phi();
2504     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2505     //fill the fDirectCorrelations:     
2506     
2507     // 2-p
2508     fDirectCorrelations->Fill(101.,cos(n*(phi2-phi3)),wPhi1*wPhi2*wPhi3); // <w1 w2 w3 cos(n(phi2-phi3))>
2509     fDirectCorrelations->Fill(102.,cos(n*(phi1-phi3)),pow(wPhi2,2.)*wPhi3); // <w2^2 w3 cos(n(psi1-phi2))>
2510     
2511     // 3-p            
2512     fDirectCorrelations->Fill(110.,cos(n*(2.*phi1-phi2-phi3)),wPhi1*wPhi2*wPhi3); // <w1 w2 w3 cos(n(2psi1-phi2-phi3))>
2513     fDirectCorrelations->Fill(111.,cos(n*(phi1+phi2-2.*phi3)),wPhi2*pow(wPhi3,2.)); // <w2 w3^2 cos(n(psi1+phi2-2.*phi3))>
2514     
2515     
2516     //fDirectCorrelations->Fill(46.,cos(n*(phi1+phi2-2.*phi3)),1);//<3'>_{n,n|2n}    
2517    }//end of for(Int_t i3=0;i3<nPrim;i3++)  
2518   }//end of for(Int_t i2=0;i2<nPrim;i2++)  
2519  }//end of for(Int_t i1=0;i1<nPrim;i1++) 
2520  
2521  
2522  */
2523  
2524  
2525  /*
2526  
2527  //<4'>_{n,n|n,n}
2528  for(Int_t i1=0;i1<nPrim;i1++)
2529  {
2530   fTrack=anEvent->GetTrack(i1);
2531   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2532   tempCounter++;
2533   phi1=fTrack->Phi();
2534   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2535   for(Int_t i2=0;i2<nPrim;i2++)
2536   {
2537    if(i2==i1)continue;
2538    fTrack=anEvent->GetTrack(i2);
2539    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2540    phi2=fTrack->Phi();
2541    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2542    for(Int_t i3=0;i3<nPrim;i3++)
2543    { 
2544     if(i3==i1||i3==i2)continue;
2545     fTrack=anEvent->GetTrack(i3);
2546     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2547     phi3=fTrack->Phi();
2548     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2549     for(Int_t i4=0;i4<nPrim;i4++)
2550     {
2551      if(i4==i1||i4==i2||i4==i3)continue;
2552      fTrack=anEvent->GetTrack(i4);
2553      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2554      phi4=fTrack->Phi();
2555      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2556      //fill the fDirectCorrelations:
2557      // 4-p            
2558      fDirectCorrelations->Fill(120.,cos(n*(phi1+phi2-phi3-phi4)),wPhi2*wPhi3*wPhi4); // <w1 w2 w3 cos(n(psi1+phi1-phi2-phi3))> 
2559      
2560        tempLoop+=wPhi2*wPhi3*wPhi4;
2561      
2562     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2563    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2564   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2565  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2566   
2567  */
2568  
2569  
2570  
2571  
2572  
2573  
2574  
2575  
2576  
2577  
2578  
2579  
2580  
2581  
2582  
2583  
2584  /*
2585  
2586  
2587  
2588  
2589  
2590  
2591  
2592  //<5'>_{2n,n|n,n,n}
2593  for(Int_t i1=0;i1<nPrim;i1++)
2594  {
2595   fTrack=anEvent->GetTrack(i1);
2596   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2597   phi1=fTrack->Phi();
2598   for(Int_t i2=0;i2<nPrim;i2++)
2599   {
2600    if(i2==i1)continue;
2601    fTrack=anEvent->GetTrack(i2);
2602    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2603    phi2=fTrack->Phi();
2604    for(Int_t i3=0;i3<nPrim;i3++)
2605    { 
2606     if(i3==i1||i3==i2)continue;
2607     fTrack=anEvent->GetTrack(i3);
2608     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2609     phi3=fTrack->Phi();
2610     for(Int_t i4=0;i4<nPrim;i4++)
2611     {
2612      if(i4==i1||i4==i2||i4==i3)continue;
2613      fTrack=anEvent->GetTrack(i4);
2614      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2615      phi4=fTrack->Phi();//
2616      for(Int_t i5=0;i5<nPrim;i5++)
2617      {
2618       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2619       fTrack=anEvent->GetTrack(i5);
2620       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2621       phi5=fTrack->Phi();    
2622       //fill the fDirectCorrelations:if(bNestedLoops)
2623       fDirectCorrelations->Fill(55.,cos(2.*n*phi1+n*phi2-n*phi3-n*phi4-n*phi5),1);//<5'>_{2n,n|n,n,n}
2624      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2625     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2626    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2627   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2628  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2629  
2630  //<6'>_{n,n,n|n,n,n}
2631  for(Int_t i1=0;i1<nPrim;i1++)
2632  {
2633   fTrack=anEvent->GetTrack(i1);
2634   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2635   phi1=fTrack->Phi();
2636   for(Int_t i2=0;i2<nPrim;i2++)
2637   {
2638    if(i2==i1)continue;
2639    fTrack=anEvent->GetTrack(i2);
2640    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2641    phi2=fTrack->Phi();
2642    for(Int_t i3=0;i3<nPrim;i3++)
2643    { 
2644     if(i3==i1||i3==i2)continue;
2645     fTrack=anEvent->GetTrack(i3);
2646     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2647     phi3=fTrack->Phi();
2648     for(Int_t i4=0;i4<nPrim;i4++)
2649     {
2650      if(i4==i1||i4==i2||i4==i3)continue;
2651      fTrack=anEvent->GetTrack(i4);
2652      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2653      phi4=fTrack->Phi();
2654      for(Int_t i5=0;i5<nPrim;i5++)
2655      {
2656       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2657       fTrack=anEvent->GetTrack(i5);
2658       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2659       phi5=fTrack->Phi();    
2660       for(Int_t i6=0;i6<nPrim;i6++)
2661       {
2662        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2663        fTrack=anEvent->GetTrack(i6);
2664        if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2665        phi6=fTrack->Phi();  
2666        //fill the fDirectCorrelations:
2667        fDirectCorrelations->Fill(60.,cos(n*(phi1+phi2+phi3-phi4-phi5-phi6)),1);//<6'>_{n,n,n|n,n,n}
2668       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
2669      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2670     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2671    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2672   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2673  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2674  
2675  //<7'>_{2n,n,n|n,n,n,n}
2676  for(Int_t i1=0;i1<nPrim;i1++)
2677  {
2678   fTrack=anEvent->GetTrack(i1);
2679   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2680   phi1=fTrack->Phi();
2681   for(Int_t i2=0;i2<nPrim;i2++)
2682   {
2683    if(i2==i1)continue;
2684    fTrack=anEvent->GetTrack(i2);
2685    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2686    phi2=fTrack->Phi();
2687    for(Int_t i3=0;i3<nPrim;i3++)
2688    { 
2689     if(i3==i1||i3==i2)continue;
2690     fTrack=anEvent->GetTrack(i3);
2691     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2692     phi3=fTrack->Phi();
2693     for(Int_t i4=0;i4<nPrim;i4++)
2694     {
2695      if(i4==i1||i4==i2||i4==i3)continue;
2696      fTrack=anEvent->GetTrack(i4);
2697      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2698      phi4=fTrack->Phi();
2699      for(Int_t i5=0;i5<nPrim;i5++)
2700      {
2701       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2702       fTrack=anEvent->GetTrack(i5);
2703       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2704       phi5=fTrack->Phi();    
2705       for(Int_t i6=0;i6<nPrim;i6++)
2706       {
2707        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2708        fTrack=anEvent->GetTrack(i6);
2709        if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2710        phi6=fTrack->Phi();
2711        for(Int_t i7=0;i7<nPrim;i7++)
2712        {
2713         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2714         fTrack=anEvent->GetTrack(i7);
2715         if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2716         phi7=fTrack->Phi();   
2717         //fill the fDirectCorrelations:
2718         fDirectCorrelations->Fill(65.,cos(2.*n*phi1+n*phi2+n*phi3-n*phi4-n*phi5-n*phi6-n*phi7),1);//<7'>_{2n,n,n|n,n,n,n}
2719        }//end of for(Int_t i7=0;i7<nPrim;i7++)  
2720       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
2721      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2722     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2723    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2724   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2725  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2726  
2727  //<8'>_{n,n,n,n|n,n,n,n}
2728  for(Int_t i1=0;i1<nPrim;i1++)
2729  {
2730   fTrack=anEvent->GetTrack(i1);
2731   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2732   phi1=fTrack->Phi();
2733   for(Int_t i2=0;i2<nPrim;i2++)
2734   {
2735    if(i2==i1)continue;
2736    fTrack=anEvent->GetTrack(i2);
2737    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2738    phi2=fTrack->Phi();
2739    for(Int_t i3=0;i3<nPrim;i3++)
2740    { 
2741     if(i3==i1||i3==i2)continue;
2742     fTrack=anEvent->GetTrack(i3);
2743     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2744     phi3=fTrack->Phi();
2745     for(Int_t i4=0;i4<nPrim;i4++)
2746     {
2747      if(i4==i1||i4==i2||i4==i3)continue;
2748      fTrack=anEvent->GetTrack(i4);
2749      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2750      phi4=fTrack->Phi();
2751      for(Int_t i5=0;i5<nPrim;i5++)
2752      {
2753       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2754       fTrack=anEvent->GetTrack(i5);
2755       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2756       phi5=fTrack->Phi();    
2757       for(Int_t i6=0;i6<nPrim;i6++)
2758       {
2759        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2760        fTrack=anEvent->GetTrack(i6);
2761        if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2762        phi6=fTrack->Phi();
2763        for(Int_t i7=0;i7<nPrim;i7++)
2764        {
2765         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2766         fTrack=anEvent->GetTrack(i7);
2767         if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2768         phi7=fTrack->Phi();
2769         for(Int_t i8=0;i8<nPrim;i8++)
2770         {
2771          if(i8==i1||i8==i2||i8==i3||i8==i4||i8==i5||i8==i6||i8==i7)continue;
2772          fTrack=anEvent->GetTrack(i8);
2773          if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2774          phi8=fTrack->Phi();           
2775          //fill the fDirectCorrelations:
2776          fDirectCorrelations->Fill(70.,cos(n*(phi1+phi2+phi3+phi4-phi5-phi6-phi7-phi8)),1);//<8'>_{n,n,n,n|n,n,n,n}
2777         }//end of for(Int_t i8=0;i8<nPrim;i8++) 
2778        }//end of for(Int_t i7=0;i7<nPrim;i7++)  
2779       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
2780      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2781     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2782    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2783   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2784  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2785  
2786  
2787  
2788  
2789  
2790  
2791  
2792  
2793  
2794  
2795  */
2796  
2797  
2798  
2799  
2800  
2801  
2802  
2803  
2804  
2805  
2806  
2807  
2808  
2809  
2810  
2811  
2812  
2813  
2814  
2815  
2816  
2817  
2818
2819  
2820  
2821  
2822  
2823  
2824  
2825  
2826  
2827  
2828 //--------------------------------------------------------------------------------------------------------------------------------  
2829
2830  //}//end of if(nPrim>0&&nPrim<12)
2831 }//end of Make()
2832
2833 //================================================================================================================
2834
2835 void AliFlowAnalysisWithQCumulants::Finish()
2836 {
2837  //calculate the final results
2838  
2839  // harmonics
2840  Int_t n = 2; // to be improved 
2841  
2842  //--------------------------------------------------------------------------------------------------------- 
2843  // avarage multiplicity
2844  Double_t AvMPOI = (fCommonHists2nd->GetHistMultDiff())->GetMean(); // to be improved 
2845  Double_t AvMRP  = (fCommonHists2nd->GetHistMultInt())->GetMean(); // to be improved 
2846
2847  // number of events
2848  Double_t nEvtsPOI = (fCommonHists2nd->GetHistMultDiff())->GetEntries(); // to be improved
2849  Double_t nEvtsRP  = (fCommonHists2nd->GetHistMultInt())->GetEntries(); // to be improved
2850  //---------------------------------------------------------------------------------------------------------
2851  
2852  //---------------------------------------------------------------------------------------------------------
2853  // 2-, 4-, 6- and 8-particle azimuthal correlation:
2854  Double_t two   = fQCorrelations->GetBinContent(1);  //<<2>>_{n|n}
2855  Double_t four  = fQCorrelations->GetBinContent(11); //<<4>>_{n,n|n,n}
2856  Double_t six   = fQCorrelations->GetBinContent(24); //<<6>>_{n,n,n|n,n,n}
2857  Double_t eight = fQCorrelations->GetBinContent(31); //<<8>>_{n,n,n,n|n,n,n,n}
2858  
2859  // 2nd, 4th, 6th and 8th order Q-cumulant:
2860  Double_t secondOrderQCumulant = two; //c_n{2} 
2861  Double_t fourthOrderQCumulant = four-2.*pow(two,2.); //c_n{4}
2862  Double_t sixthOrderQCumulant  = six-9.*two*four+12.*pow(two,3.); //c_n{6}
2863  Double_t eightOrderQCumulant  = eight-16.*two*six-18.*pow(four,2.)+144.*pow(two,2.)*four-144.*pow(two,4.); //c_n{8} 
2864  
2865  // "no-name" integrated flow estimates from Q-cumulants:
2866  Double_t vn2=0.,vn4=0.,vn6=0.,vn8=0.;
2867  // Double_t sd2=0.,sd4=0.,sd6=0.,sd8=0.; 
2868  Double_t sd6=0.,sd8=0.; // to be improved/removed
2869  if(secondOrderQCumulant>0.)
2870  {
2871   vn2 = pow(secondOrderQCumulant,0.5); //v_n{2}
2872  } 
2873  if(fourthOrderQCumulant<0.)
2874  {
2875   vn4 = pow(-fourthOrderQCumulant,1./4.); //v_n{4}
2876  } 
2877  if(sixthOrderQCumulant>0.)
2878  {
2879   vn6 = pow((1./4.)*sixthOrderQCumulant,1./6.); //v_n{6}
2880  } 
2881  if(eightOrderQCumulant<0.)
2882  {
2883   vn8 = pow((-1./33.)*eightOrderQCumulant,1./8.); //v_n{8}
2884  }
2885  //---------------------------------------------------------------------------------------------------------
2886
2887  //--------------------------------------------------------------------------------------------------------- 
2888  // weighted 2-, 4-, 6- and 8-particle azimuthal correlation:
2889  Double_t twoW   = fWeightedQCorrelations->GetBinContent(1);  //<<2>>_{n|n}
2890  Double_t fourW  = fWeightedQCorrelations->GetBinContent(21); //<<4>>_{n,n|n,n}
2891  // Double_t sixW   = fWeightedQCorrelations->GetBinContent(24); //<<6>>_{n,n,n|n,n,n}
2892  // Double_t eightW = fWeightedQCorrelations->GetBinContent(31); //<<8>>_{n,n,n,n|n,n,n,n}
2893  
2894  // 2nd, 4th, 6th and 8th order weighted Q-cumulant:
2895  Double_t secondOrderQCumulantW = twoW; //c_n{2} 
2896  Double_t fourthOrderQCumulantW = fourW-2.*pow(twoW,2.); //c_n{4}
2897  // Double_t sixthOrderQCumulantW  = sixW-9.*twoW*fourW+12.*pow(twoW,3.); //c_n{6}
2898  // Double_t eightOrderQCumulantW  = eightW-16.*twoW*sixW-18.*pow(fourW,2.)+144.*pow(twoW,2.)*fourW-144.*pow(twoW,4.); //c_n{8} 
2899  
2900  // "no-name" integrated flow estimates from weighted Q-cumulants:
2901  cout<<endl;
2902  cout<<"**************************************"<<endl;
2903  cout<<"**************************************"<<endl;
2904  cout<<"flow estimates from Q-cumulants :"<<endl;
2905  cout<<endl;
2906  
2907  Double_t vn2W=0.,vn4W=0.;
2908  Double_t sd2W=0.,sd4W=0.; 
2909  if(secondOrderQCumulantW>0.){
2910   vn2W = pow(secondOrderQCumulantW,0.5); // weighted v_n{2}
2911   //sd2W = 0.5*pow(secondOrderQCumulantW,-0.5)*secondOrderQCumulantErrorW; // to be improved (correct treatment of errors needed)
2912   cout<<" v_"<<n<<"{2} = "<<vn2W<<" +/- "<<sd2W<<endl;
2913   fIntFlowResultsQC->SetBinContent(1,vn2W);
2914   //fIntFlowResultsQC->SetBinError(1,sd2W);
2915   //common histograms:
2916   fCommonHistsResults2nd->FillIntegratedFlow(vn2W,sd2W);
2917   //fCommonHistsResults2nd->FillChi(vn2W*pow(AvM,0.5)); // to be removed
2918  }else{
2919   cout<<" v_"<<n<<"{2} = Im"<<endl;
2920  }          
2921  if(fourW!=0. && fourthOrderQCumulantW<0.){
2922   vn4W = pow(-fourthOrderQCumulantW,1./4.); //v_n{4}
2923   //sd4W = 0.25*pow(-fourthOrderQCumulantW,-3./4.)*fourthOrderQCumulantErrorW; // to be improved (correct treatment of errors needed)
2924   cout<<" v_"<<n<<"{4} = "<<vn4W<<" +/- "<<sd4W<<endl;
2925   fIntFlowResultsQC->SetBinContent(2,vn4W);
2926   //fIntFlowResultsQC->SetBinError(2,sd4W);
2927   //common histograms:
2928   fCommonHistsResults4th->FillIntegratedFlow(vn4W,sd4W);
2929   //fCommonHistsResults4th->FillChi(vn4W*pow(AvM,0.5)); // to be removed
2930  }else{
2931   cout<<" v_"<<n<<"{4} = Im"<<endl;
2932  }
2933  // !!! to be improved (6th and 8th order are without weights) !!!
2934  if(six!=0. && sixthOrderQCumulant>0.){
2935   vn6 = pow((1./4.)*sixthOrderQCumulant,1./6.); //v_n{6}
2936   //sd6 = (1./6.)*pow(2.,-1./3.)*pow(sixthOrderQCumulant,-5./6.)*sixthOrderQCumulantError;
2937   cout<<" v_"<<n<<"{6} = "<<vn6<<" +/- "<<sd6<<endl;
2938   fIntFlowResultsQC->SetBinContent(3,vn6);
2939   //fIntFlowResultsQC->SetBinError(3,sd6);
2940   //common histograms:
2941   fCommonHistsResults6th->FillIntegratedFlow(vn6,sd6);
2942   //fCommonHistsResults6th->FillChi(vn6*pow(AvM,0.5));//to be removed
2943  }else{
2944   cout<<" v_"<<n<<"{6} = Im"<<endl;
2945  }
2946  if(eight!=0. && eightOrderQCumulant<0.){
2947   vn8 = pow((-1./33.)*eightOrderQCumulant,1./8.); //v_n{8}
2948   cout<<" v_"<<n<<"{8} = "<<vn8<<" +/- "<<sd8<<endl;
2949   fIntFlowResultsQC->SetBinContent(4,vn8);
2950   //fIntFlowResultsQC->SetBinError(4,sd8);
2951   //common histograms:
2952   fCommonHistsResults8th->FillIntegratedFlow(vn8,sd8);
2953   //fCommonHistsResults8th->FillChi(vn8*pow(AvM,0.5));//to be removed
2954  }else{
2955   cout<<" v_"<<n<<"{8} = Im"<<endl;
2956  }
2957  cout<<endl;
2958  cout<<"   nEvts = "<<nEvtsRP<<", AvM = "<<AvMRP<<endl; // to be improved
2959  cout<<"**************************************"<<endl;
2960  cout<<"**************************************"<<endl;
2961  cout<<endl; 
2962 //--------------------------------------------------------------------------------------------------------- 
2963  
2964 //---------------------------------------------------------------------------------------------------------
2965 // differential flow (POI)
2966 Int_t nBinsPtPOI = f2WPerPtBin1n1nPOI->GetNbinsX();
2967 Int_t nBinsEtaPOI = f4WPerPtBin1n1n1n1nPOI->GetNbinsX();
2968
2969 // Pt:
2970 Double_t secondOrderQCumulantDiffFlowPtPOIW = 0.;
2971 Double_t fourthOrderQCumulantDiffFlowPtPOIW = 0.;
2972
2973 Double_t dVn2ndPOIW=0.,dSd2ndPOIW=0.,dDiffvn2ndPOIW=0.,dYield2ndPOIW=0.,dSum2ndPOIW=0.;
2974 Double_t dVn4thPOIW=0.,dSd4thPOIW=0.,dDiffvn4thPOIW=0.,dYield4thPOIW=0.,dSum4thPOIW=0.;
2975
2976 for(Int_t bb=1;bb<nBinsPtPOI+1;bb++)
2977 {
2978  // QC{2}
2979  if(f2WPerPtBin1n1nPOI->GetBinEntries(bb)>0.&&vn2W!=0)
2980  { 
2981   secondOrderQCumulantDiffFlowPtPOIW = f2WPerPtBin1n1nPOI->GetBinContent(bb); // with weights
2982   fDiffFlowResults2ndOrderQC->SetBinContent(bb,secondOrderQCumulantDiffFlowPtPOIW/vn2W);
2983   // common histogram:
2984   fCommonHistsResults2nd->FillDifferentialFlowPtPOI(bb,secondOrderQCumulantDiffFlowPtPOIW/vn2W, 0.); //to be improved (errors && bb or bb+1 ?)
2985   // -------------------------------------------------------------------
2986   // integrated flow (weighted, POI, Pt, 2nd order):
2987   dDiffvn2ndPOIW=(fCommonHistsResults2nd->GetHistDiffFlowPtPOI())->GetBinContent(bb);
2988   dYield2ndPOIW=(fCommonHists2nd->GetHistPtDiff())->GetBinContent(bb);
2989   dVn2ndPOIW+=dDiffvn2ndPOIW*dYield2ndPOIW;
2990   dSum2ndPOIW+=dYield2ndPOIW;
2991   // -------------------------------------------------------------------
2992  }
2993  // QC{4]
2994  if(f4WPerPtBin1n1n1n1nPOI->GetBinEntries(bb)>0.&&vn4W!=0.)
2995  {
2996   fourthOrderQCumulantDiffFlowPtPOIW = f4WPerPtBin1n1n1n1nPOI->GetBinContent(bb)-2.*f2WPerPtBin1n1nPOI->GetBinContent(bb)*pow(vn2W,2.); // with weights
2997   fDiffFlowResults4thOrderQC->SetBinContent(bb,-1.*fourthOrderQCumulantDiffFlowPtPOIW/pow(vn4W,3.));
2998   //common histogram:
2999   fCommonHistsResults4th->FillDifferentialFlowPtPOI(bb,-1.*fourthOrderQCumulantDiffFlowPtPOIW/pow(vn4W,3.), 0.); //to be improved (errors)
3000   // -------------------------------------------------------------------
3001   //integrated flow (POI, Pt, 4th order):
3002   dDiffvn4thPOIW=(fCommonHistsResults4th->GetHistDiffFlowPtPOI())->GetBinContent(bb);
3003   dYield4thPOIW=(fCommonHists4th->GetHistPtDiff())->GetBinContent(bb);
3004   dVn4thPOIW+=dDiffvn4thPOIW*dYield4thPOIW;
3005   dSum4thPOIW+=dYield4thPOIW;
3006   // -------------------------------------------------------------------
3007  }
3008 }      
3009
3010 cout<<endl;
3011 cout<<"**************************************"<<endl;
3012 cout<<"**************************************"<<endl;
3013 cout<<"flow estimates from Q-cumulants (POI):"<<endl;
3014 cout<<endl;
3015 //storing the final results for integrated flow (POI):
3016 // QC{2}
3017 if(dSum2ndPOIW && fCommonHistsResults2nd)
3018 {
3019  dVn2ndPOIW/=dSum2ndPOIW;
3020  fCommonHistsResults2nd->FillIntegratedFlowPOI(dVn2ndPOIW,0.); // to be improved (errors)
3021  cout<<" v_"<<n<<"{2} = "<<dVn2ndPOIW<<" +/- "<<dSd2ndPOIW<<endl;
3022 }else 
3023  {
3024   cout<<" v_"<<n<<"{2} = Im"<<endl;
3025  }
3026
3027 // QC{4}
3028 if(dSum4thPOIW && fCommonHistsResults4th)
3029 {
3030  dVn4thPOIW/=dSum4thPOIW;
3031  fCommonHistsResults4th->FillIntegratedFlowPOI(dVn4thPOIW,0.); // to be improved (errors)
3032  cout<<" v_"<<n<<"{4} = "<<dVn4thPOIW<<" +/- "<<dSd4thPOIW<<endl;
3033 }else
3034  {
3035   cout<<" v_"<<n<<"{4} = Im"<<endl;
3036  }
3037
3038 cout<<endl;
3039 cout<<"   nEvts = "<<nEvtsPOI<<", AvM = "<<AvMPOI<<endl;
3040 cout<<"**************************************"<<endl;
3041 cout<<"**************************************"<<endl;
3042 cout<<endl;  
3043
3044 //Eta:
3045 Double_t secondOrderQCumulantDiffFlowEtaPOIW = 0.;
3046 Double_t fourthOrderQCumulantDiffFlowEtaPOIW = 0.;
3047
3048 for(Int_t bb=1;bb<nBinsEtaPOI+1;bb++)
3049 {
3050  if(f2WPerEtaBin1n1nPOI->GetBinEntries(bb)>0.&&vn2W!=0)
3051  {
3052   secondOrderQCumulantDiffFlowEtaPOIW = f2WPerEtaBin1n1nPOI->GetBinContent(bb); // with weights
3053   fDiffFlowResults2ndOrderQC->SetBinContent(bb,secondOrderQCumulantDiffFlowEtaPOIW/vn2W);
3054   //common histogram:
3055   fCommonHistsResults2nd->FillDifferentialFlowEtaPOI(bb,secondOrderQCumulantDiffFlowEtaPOIW/vn2W, 0.);//to be improved (errors)
3056  }
3057  if(f4WPerEtaBin1n1n1n1nPOI->GetBinEntries(bb)>0.&&vn4W!=0.)
3058  {
3059   fourthOrderQCumulantDiffFlowEtaPOIW = f4WPerEtaBin1n1n1n1nPOI->GetBinContent(bb)-2.*f2WPerEtaBin1n1nPOI->GetBinContent(bb)*pow(vn2W,2.); // with weights
3060   fDiffFlowResults4thOrderQC->SetBinContent(bb,-1.*fourthOrderQCumulantDiffFlowEtaPOIW/pow(vn4W,3.));
3061   //common histogram:
3062   fCommonHistsResults4th->FillDifferentialFlowEtaPOI(bb,-1.*fourthOrderQCumulantDiffFlowEtaPOIW/pow(vn4W,3.), 0.);//to be improved (errors)
3063  }
3064 }      
3065 //------------------------------------------------------------
3066  
3067  
3068  
3069  
3070  
3071  //---------------------------------------------------------------------------------------------------------
3072 // differential flow (RP)
3073 Int_t nBinsPtRP = f2WPerPtBin1n1nRP->GetNbinsX();
3074 Int_t nBinsEtaRP = f4WPerPtBin1n1n1n1nRP->GetNbinsX();
3075
3076 // Pt:
3077 Double_t secondOrderQCumulantDiffFlowPtRPW = 0.;
3078 Double_t fourthOrderQCumulantDiffFlowPtRPW = 0.;
3079
3080 Double_t dVn2ndRPW=0.,dSd2ndRPW=0.,dDiffvn2ndRPW=0.,dYield2ndRPW=0.,dSum2ndRPW=0.;
3081 Double_t dVn4thRPW=0.,dSd4thRPW=0.,dDiffvn4thRPW=0.,dYield4thRPW=0.,dSum4thRPW=0.;
3082
3083 for(Int_t bb=1;bb<nBinsPtRP+1;bb++)
3084 {
3085  // QC{2}
3086  if(f2WPerPtBin1n1nRP->GetBinEntries(bb)>0.&&vn2W!=0)
3087  { 
3088   secondOrderQCumulantDiffFlowPtRPW = f2WPerPtBin1n1nRP->GetBinContent(bb); // with weights
3089   fDiffFlowResults2ndOrderQC->SetBinContent(bb,secondOrderQCumulantDiffFlowPtRPW/vn2W);
3090   // common histogram:
3091   fCommonHistsResults2nd->FillDifferentialFlowPtRP(bb,secondOrderQCumulantDiffFlowPtRPW/vn2W, 0.); //to be improved (errors && bb or bb+1 ?)
3092   // -------------------------------------------------------------------
3093   // integrated flow (weighted, RP, Pt, 2nd order):
3094   dDiffvn2ndRPW=(fCommonHistsResults2nd->GetHistDiffFlowPtRP())->GetBinContent(bb);
3095   dYield2ndRPW=(fCommonHists2nd->GetHistPtInt())->GetBinContent(bb);
3096   dVn2ndRPW+=dDiffvn2ndRPW*dYield2ndRPW;
3097   dSum2ndRPW+=dYield2ndRPW;
3098   // -------------------------------------------------------------------
3099  }
3100  // QC{4]
3101  if(f4WPerPtBin1n1n1n1nRP->GetBinEntries(bb)>0.&&vn4W!=0.)
3102  {
3103   fourthOrderQCumulantDiffFlowPtRPW = f4WPerPtBin1n1n1n1nRP->GetBinContent(bb)-2.*f2WPerPtBin1n1nRP->GetBinContent(bb)*pow(vn2W,2.); // with weights
3104   fDiffFlowResults4thOrderQC->SetBinContent(bb,-1.*fourthOrderQCumulantDiffFlowPtRPW/pow(vn4W,3.));
3105   //common histogram:
3106   fCommonHistsResults4th->FillDifferentialFlowPtRP(bb,-1.*fourthOrderQCumulantDiffFlowPtRPW/pow(vn4W,3.), 0.); //to be improved (errors)
3107   // -------------------------------------------------------------------
3108   //integrated flow (RP, Pt, 4th order):
3109   dDiffvn4thRPW=(fCommonHistsResults4th->GetHistDiffFlowPtRP())->GetBinContent(bb);
3110   dYield4thRPW=(fCommonHists4th->GetHistPtInt())->GetBinContent(bb);
3111   dVn4thRPW+=dDiffvn4thRPW*dYield4thRPW;
3112   dSum4thRPW+=dYield4thRPW;
3113   // -------------------------------------------------------------------
3114  }
3115 }      
3116
3117 cout<<endl;
3118 cout<<"**************************************"<<endl;
3119 cout<<"**************************************"<<endl;
3120 cout<<"flow estimates from Q-cumulants (RP):"<<endl;
3121 cout<<endl;
3122 //storing the final results for integrated flow (RP):
3123 // QC{2}
3124 if(dSum2ndRPW && fCommonHistsResults2nd)
3125 {
3126  dVn2ndRPW/=dSum2ndRPW;
3127  fCommonHistsResults2nd->FillIntegratedFlowRP(dVn2ndRPW,0.); // to be improved (errors)
3128  cout<<" v_"<<n<<"{2} = "<<dVn2ndRPW<<" +/- "<<dSd2ndRPW<<endl;
3129 }else 
3130  {
3131   cout<<" v_"<<n<<"{2} = Im"<<endl;
3132  }
3133
3134 // QC{4}
3135 if(dSum4thRPW && fCommonHistsResults4th)
3136 {
3137  dVn4thRPW/=dSum4thRPW;
3138  fCommonHistsResults4th->FillIntegratedFlowRP(dVn4thRPW,0.); // to be improved (errors)
3139  cout<<" v_"<<n<<"{4} = "<<dVn4thRPW<<" +/- "<<dSd4thRPW<<endl;
3140 }else
3141  {
3142   cout<<" v_"<<n<<"{4} = Im"<<endl;
3143  }
3144
3145 cout<<endl;
3146 cout<<"   nEvts = "<<nEvtsRP<<", AvM = "<<AvMRP<<endl;
3147 cout<<"**************************************"<<endl;
3148 cout<<"**************************************"<<endl;
3149 cout<<endl;  
3150
3151 //Eta:
3152 Double_t secondOrderQCumulantDiffFlowEtaRPW = 0.;
3153 Double_t fourthOrderQCumulantDiffFlowEtaRPW = 0.;
3154
3155 for(Int_t bb=1;bb<nBinsEtaRP+1;bb++)
3156 {
3157  if(f2WPerEtaBin1n1nRP->GetBinEntries(bb)>0.&&vn2W!=0)
3158  {
3159   secondOrderQCumulantDiffFlowEtaRPW = f2WPerEtaBin1n1nRP->GetBinContent(bb); // with weights
3160   fDiffFlowResults2ndOrderQC->SetBinContent(bb,secondOrderQCumulantDiffFlowEtaRPW/vn2W);
3161   //common histogram:
3162   fCommonHistsResults2nd->FillDifferentialFlowEtaRP(bb,secondOrderQCumulantDiffFlowEtaRPW/vn2W, 0.);//to be improved (errors)
3163  }
3164  if(f4WPerEtaBin1n1n1n1nRP->GetBinEntries(bb)>0.&&vn4W!=0.)
3165  {
3166   fourthOrderQCumulantDiffFlowEtaRPW = f4WPerEtaBin1n1n1n1nRP->GetBinContent(bb)-2.*f2WPerEtaBin1n1nRP->GetBinContent(bb)*pow(vn2W,2.); // with weights
3167   fDiffFlowResults4thOrderQC->SetBinContent(bb,-1.*fourthOrderQCumulantDiffFlowEtaRPW/pow(vn4W,3.));
3168   //common histogram:
3169   fCommonHistsResults4th->FillDifferentialFlowEtaRP(bb,-1.*fourthOrderQCumulantDiffFlowEtaRPW/pow(vn4W,3.), 0.);//to be improved (errors)
3170  }
3171 }      
3172 //------------------------------------------------------------
3173  
3174  
3175  
3176  /*
3177  cout<<"0.5 < Pt < 0.6 GeV"<<endl;                                
3178  cout<<endl;                                       
3179  cout<<"<w2 cos(n(psi1-phi2))> from Q-vectors                    = "<<f2WPerPtBin1n1nPOI->GetBinContent(6)<<endl;
3180  //cout<<"<w2 cos(n(psi1-phi2))> from Q-vectors                    = "<<f2WPerPtBin1n1nRP->GetName()<<endl;
3181  cout<<"<w2 cos(n(psi1-phi2))> from Q-vectors                    = "<<fDirectCorrelations->GetBinContent(101)<<endl;
3182  cout<<endl;  
3183  cout<<"<w2 w3 w4 cos(n(psi1+phi2-phi3-phi4))> from Q-vectors    = "<<f4WPerPtBin1n1n1n1nPOI->GetBinContent(6)<<endl;
3184  //cout<<"<w2 w3 w4 cos(n(psi1+phi2-phi3-phi4))> from Q-vectors    = "<<f4WPerPtBin1n1n1n1nRP->GetName()<<endl;
3185  cout<<"<w2 w3 w4 cos(n(psi1+phi2-phi3-phi4))> from nested loops = "<<fDirectCorrelations->GetBinContent(121)<<endl;   
3186  cout<<endl; 
3187  
3188  
3189  
3190  
3191  
3192   */
3193  
3194  
3195  
3196  /*
3197    
3198     
3199  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
3200  //                   !!!! to be removed !!!!  
3201  
3202  TCanvas* qvectorPlot = new TCanvas("qvectorPlot","Q-vector Plot",1000,1000);
3203  
3204  qvectorPlot->cd(1);
3205  
3206  TH1D* style = new TH1D("style","Q-vectors",100,-244,244); 
3207  (style->GetYaxis())->SetRangeUser(-244,244);
3208  
3209  style->Draw();
3210   
3211  Int_t nBins=fQvectorForEachEventX->GetNbinsX();
3212  Double_t qxxx=0.,qyyy=0.;
3213  //cout<<"nBins = "<<nBins<<endl;   
3214  //cout<<fQvectorForEachEventX->GetBinEntries(4)<<endl;    
3215  //cout<<fQvectorForEachEventY->GetBinEntries(4)<<endl;     
3216  
3217  for(Int_t b=1;b<nBins+1;b++)
3218  {
3219   if(fQvectorForEachEventX->GetBinEntries(b)==1 && fQvectorForEachEventY->GetBinEntries(b)==1)
3220   {
3221    qxxx=fQvectorForEachEventX->GetBinContent(b);
3222    qyyy=fQvectorForEachEventY->GetBinContent(b);
3223    //cout<<qxxx<<" "<<qyyy<<endl;
3224    TArrow *qvector = new TArrow(0.0,0.0,qxxx,qyyy,0.0144,"|>");
3225    qvector->SetAngle(40);
3226    qvector->SetLineWidth(2);
3227    qvector->Draw("");
3228   }
3229  }  
3230  
3231  //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx       
3232
3233   */     
3234                            
3235 }
3236
3237 //================================================================================================================
3238
3239 void AliFlowAnalysisWithQCumulants::WriteHistograms(TString outputFileName)
3240 {
3241  //store the final results in output .root file
3242  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
3243  output->WriteObject(fHistList, "cobjQC","SingleKey");
3244  delete output;
3245 }
3246
3247 //================================================================================================================
3248
3249
3250