forgotten delete
[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  delete ptReq1nPrimePrime;
1931  delete ptImq1nPrimePrime;
1932  delete ptReq2nPrimePrime; 
1933  delete ptImq2nPrimePrime;
1934
1935  delete req1nW2PrimePrimePt;
1936  delete imq1nW2PrimePrimePt;
1937  delete req2nW1PrimePrimePt;
1938  delete imq2nW1PrimePrimePt;
1939  
1940  delete sumOfW1upTomPrimePrimePt;
1941  delete sumOfW2upTomPrimePrimePt;
1942  delete sumOfW3upTomPrimePrimePt;
1943  
1944  delete etaReq1nPrimePrime;
1945  delete etaImq1nPrimePrime;
1946  delete etaReq2nPrimePrime; 
1947  delete etaImq2nPrimePrime;
1948  
1949  delete req1nW2PrimePrimeEta;
1950  delete imq1nW2PrimePrimeEta;
1951  delete req2nW1PrimePrimeEta;
1952  delete imq2nW1PrimePrimeEta;
1953  
1954  delete sumOfW1upTomPrimePrimeEta;
1955  delete sumOfW2upTomPrimePrimeEta;
1956  delete sumOfW3upTomPrimePrimeEta;
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
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106  /*
2107
2108  //if(bNestedLoops)to be improved
2109  //{ to be improved               
2110  //-------------------------------------------------------------------------------------------------------------------------------- 
2111  //
2112  //                                          **********************
2113  //                                          **** NESTED LOOPS ****
2114  //                                          ********************** 
2115  //
2116  // Remark 1: multi-particle correlations calculated with nested loops are stored in fDirectCorrelations.
2117  // 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)
2118  //
2119  // binning details of fDirectCorrelations (integrated flow):
2120  //..........................................................................
2121  // 1st bin: weighted <2>_{n|n}   = <w1   w2   cos( n*(phi1-phi2))>
2122  // 2nd bin: weighted <2>_{2n|2n} = <w1^2 w2^2 cos(2n*(phi1-phi2))>
2123  // 3rd bin: weighted <2>_{3n|3n} = <w1^3 w2^3 cos(3n*(phi1-phi2))>
2124  // 4th bin: weighted <2>_{4n|4n} = <w1^4 w2^4 cos(4n*(phi1-phi2))>
2125  // 5th bin: weighted <2>_{n|n} = <w1^3 w2 cos(n*(phi1-phi2))>
2126  
2127  // 11th bin: weighted <3>_{2n|n,n} = <w1^2 w2 w3 cos(n*(2phi1-phi2-phi3))>
2128  //..........................................................................
2129  
2130
2131  Double_t phi1=0., phi2=0., phi3=0., phi4=0., phi5=0., phi6=0., phi7=0., phi8=0.;
2132  Double_t wPhi1=1., wPhi2=1., wPhi3=1., wPhi4=1., wPhi5=1., wPhi6=1., wPhi7=1., wPhi8=1.;
2133
2134
2135
2136  
2137
2138
2139  Double_t tempLoop = 0.;
2140  Int_t tempCounter = 0;
2141  
2142  
2143  
2144  
2145  
2146  for(Int_t i1=0;i1<dMult;i1++)
2147  {
2148   fTrack=anEvent->GetTrack(i1);
2149   phi1=fTrack->Phi();
2150   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2151   for(Int_t i2=0;i2<dMult;i2++)
2152   {
2153    if(i2==i1)continue;
2154    fTrack=anEvent->GetTrack(i2);
2155    phi2=fTrack->Phi();
2156    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2157    // 2-p
2158    fDirectCorrelations->Fill(0.,cos(n*(phi1-phi2)),wPhi1*wPhi2);                  // <w1   w2   cos( n*(phi1-phi2))>
2159    fDirectCorrelations->Fill(1.,cos(2.*n*(phi1-phi2)),pow(wPhi1,2)*pow(wPhi2,2)); // <w1^2 w2^2 cos(2n*(phi1-phi2))>
2160    fDirectCorrelations->Fill(2.,cos(3.*n*(phi1-phi2)),pow(wPhi1,3)*pow(wPhi2,3)); // <w1^3 w2^3 cos(3n*(phi1-phi2))>
2161    fDirectCorrelations->Fill(3.,cos(4.*n*(phi1-phi2)),pow(wPhi1,4)*pow(wPhi2,4)); // <w1^4 w2^4 cos(4n*(phi1-phi2))> 
2162    fDirectCorrelations->Fill(4.,cos(n*(phi1-phi2)),pow(wPhi1,3)*wPhi2); // <w1^3 w2 cos(n*(phi1-phi2))>
2163   }
2164  }  
2165  
2166
2167
2168
2169
2170
2171  for(Int_t i1=0;i1<dMult;i1++)
2172  {
2173   fTrack=anEvent->GetTrack(i1);
2174   phi1=fTrack->Phi();
2175   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2176   for(Int_t i2=0;i2<dMult;i2++)
2177   {
2178    if(i2==i1)continue;
2179    fTrack=anEvent->GetTrack(i2);
2180    phi2=fTrack->Phi();
2181    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2182    for(Int_t i3=0;i3<dMult;i3++)
2183    {
2184     if(i3==i1||i3==i2)continue;
2185     fTrack=anEvent->GetTrack(i3);
2186     phi3=fTrack->Phi();
2187     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2188     // 2-p
2189     fDirectCorrelations->Fill(5.,cos(n*(phi1-phi2)),wPhi1*wPhi2*pow(wPhi3,2)); // <w1 w2 w3^2 cos(n*(phi1-phi2))>
2190     
2191     // 3-p
2192     fDirectCorrelations->Fill(10.,cos(2.*n*phi1-n*(phi2+phi3)),pow(wPhi1,2)*wPhi2*wPhi3); // <w1^2 w2 w3 cos(n*(2phi1-phi2-phi3))>
2193     
2194     
2195     fDirectCorrelations->Fill(6.,cos(3.*n*phi1-2.*n*phi2-n*phi3),pow(wPhi1,3)*pow(wPhi2,2)*wPhi3);           //<3>_{3n|2n,n}
2196     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}
2197     fDirectCorrelations->Fill(8.,cos(4.*n*phi1-3.*n*phi2-n*phi3),pow(wPhi1,4)*pow(wPhi2,3)*wPhi3);           //<3>_{4n|3n,n}
2198     
2199     
2200    }
2201   }
2202  }
2203  
2204
2205  
2206  
2207
2208
2209   
2210  //<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} 
2211  for(Int_t i1=0;i1<dMult;i1++)
2212  {
2213   fTrack=anEvent->GetTrack(i1);
2214   phi1=fTrack->Phi();
2215   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2216   for(Int_t i2=0;i2<dMult;i2++)
2217   {
2218    if(i2==i1)continue;
2219    fTrack=anEvent->GetTrack(i2);
2220    phi2=fTrack->Phi();
2221    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2222    for(Int_t i3=0;i3<dMult;i3++)
2223    {
2224     if(i3==i1||i3==i2)continue;
2225     fTrack=anEvent->GetTrack(i3);
2226     phi3=fTrack->Phi();
2227     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2228     for(Int_t i4=0;i4<dMult;i4++)
2229     {
2230      if(i4==i1||i4==i2||i4==i3)continue;
2231      fTrack=anEvent->GetTrack(i4);
2232      phi4=fTrack->Phi();
2233      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2234      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))>
2235      
2236     
2237      fDirectCorrelations->Fill(11.,cos(2.*n*phi1+n*phi2-2.*n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);      //<4>_{2n,n|2n,n}
2238      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}
2239      fDirectCorrelations->Fill(13.,cos(3.*n*phi1-n*phi2-n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);         //<4>_{3n|n,n,n}
2240      fDirectCorrelations->Fill(14.,cos(3.*n*phi1+n*phi2-3.*n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);      //<4>_{3n,n|3n,n}   
2241      fDirectCorrelations->Fill(15.,cos(3.*n*phi1+n*phi2-2.*n*phi3-2.*n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);   //<4>_{3n,n|2n,2n}
2242      fDirectCorrelations->Fill(16.,cos(4.*n*phi1-2.*n*phi2-n*phi3-n*phi4),wPhi1*wPhi2*wPhi3*wPhi4);      //<4>_{4n|2n,n,n}
2243     
2244     }  
2245    }
2246   }
2247  }
2248  
2249  
2250  
2251
2252  
2253  //<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}
2254  for(Int_t i1=0;i1<dMult;i1++)
2255  {
2256   //cout<<"i1 = "<<i1<<endl;
2257   fTrack=anEvent->GetTrack(i1);
2258   phi1=fTrack->Phi();
2259   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2260   for(Int_t i2=0;i2<dMult;i2++)
2261   {
2262    if(i2==i1)continue;
2263    fTrack=anEvent->GetTrack(i2);
2264    phi2=fTrack->Phi();
2265    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2266    for(Int_t i3=0;i3<dMult;i3++)
2267    {
2268     if(i3==i1||i3==i2)continue;
2269     fTrack=anEvent->GetTrack(i3);
2270     phi3=fTrack->Phi();
2271     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2272     for(Int_t i4=0;i4<dMult;i4++)
2273     {
2274      if(i4==i1||i4==i2||i4==i3)continue;
2275      fTrack=anEvent->GetTrack(i4);
2276      phi4=fTrack->Phi();
2277      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2278      for(Int_t i5=0;i5<dMult;i5++)
2279      {
2280       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2281       fTrack=anEvent->GetTrack(i5);
2282       phi5=fTrack->Phi();
2283       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2284       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}
2285       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}
2286       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}
2287       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}
2288      }
2289     }  
2290    }
2291   }
2292  }
2293  
2294  
2295  //<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}
2296  for(Int_t i1=0;i1<dMult;i1++)
2297  {
2298   //cout<<"i1 = "<<i1<<endl;
2299   fTrack=anEvent->GetTrack(i1);
2300   phi1=fTrack->Phi();
2301   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2302   for(Int_t i2=0;i2<dMult;i2++)
2303   {
2304    if(i2==i1)continue;
2305    fTrack=anEvent->GetTrack(i2);
2306    phi2=fTrack->Phi();
2307    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2308    for(Int_t i3=0;i3<dMult;i3++)
2309    {
2310     if(i3==i1||i3==i2)continue;
2311     fTrack=anEvent->GetTrack(i3);
2312     phi3=fTrack->Phi();
2313     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2314     for(Int_t i4=0;i4<dMult;i4++)
2315     {
2316      if(i4==i1||i4==i2||i4==i3)continue;
2317      fTrack=anEvent->GetTrack(i4);
2318      phi4=fTrack->Phi();
2319      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2320      for(Int_t i5=0;i5<dMult;i5++)
2321      {
2322       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2323       fTrack=anEvent->GetTrack(i5);
2324       phi5=fTrack->Phi();
2325       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2326       for(Int_t i6=0;i6<dMult;i6++)
2327       {
2328        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2329        fTrack=anEvent->GetTrack(i6);
2330        phi6=fTrack->Phi(); 
2331        if(phiWeights) wPhi6 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi6*nBinsPhi/TMath::TwoPi())));
2332        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}
2333        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}
2334        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}
2335        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}     
2336       } 
2337      }
2338     }  
2339    }
2340   }
2341  }
2342
2343
2344  
2345  //<7>_{2n,n,n|n,n,n,n}
2346  for(Int_t i1=0;i1<dMult;i1++)
2347  {
2348   //cout<<"i1 = "<<i1<<endl;
2349   fTrack=anEvent->GetTrack(i1);
2350   phi1=fTrack->Phi();
2351   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2352   for(Int_t i2=0;i2<dMult;i2++)
2353   {
2354    if(i2==i1)continue;
2355    fTrack=anEvent->GetTrack(i2);
2356    phi2=fTrack->Phi();
2357    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2358    for(Int_t i3=0;i3<dMult;i3++)
2359    {
2360     if(i3==i1||i3==i2)continue;
2361     fTrack=anEvent->GetTrack(i3);
2362     phi3=fTrack->Phi();
2363     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2364     for(Int_t i4=0;i4<dMult;i4++)
2365     {
2366      if(i4==i1||i4==i2||i4==i3)continue;
2367      fTrack=anEvent->GetTrack(i4);
2368      phi4=fTrack->Phi();
2369      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2370      for(Int_t i5=0;i5<dMult;i5++)
2371      {
2372       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2373       fTrack=anEvent->GetTrack(i5);
2374       phi5=fTrack->Phi();
2375       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2376       for(Int_t i6=0;i6<dMult;i6++)
2377       {
2378        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2379        fTrack=anEvent->GetTrack(i6);
2380        phi6=fTrack->Phi(); 
2381        if(phiWeights) wPhi6 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi6*nBinsPhi/TMath::TwoPi())));
2382        for(Int_t i7=0;i7<dMult;i7++)
2383        {
2384         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2385         fTrack=anEvent->GetTrack(i7);
2386         phi7=fTrack->Phi(); 
2387         if(phiWeights) wPhi7 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi7*nBinsPhi/TMath::TwoPi())));
2388         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}
2389        } 
2390       } 
2391      }
2392     }  
2393    }
2394   }
2395  }
2396  
2397  
2398  
2399  //<8>_{n,n,n,n|n,n,n,n}
2400  for(Int_t i1=0;i1<dMult;i1++)
2401  {
2402   cout<<"i1 = "<<i1<<endl;
2403   fTrack=anEvent->GetTrack(i1);
2404   phi1=fTrack->Phi();
2405   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2406   for(Int_t i2=0;i2<dMult;i2++)
2407   {
2408    if(i2==i1)continue;
2409    fTrack=anEvent->GetTrack(i2);
2410    phi2=fTrack->Phi();
2411    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2412    for(Int_t i3=0;i3<dMult;i3++)
2413    {
2414     if(i3==i1||i3==i2)continue;
2415     fTrack=anEvent->GetTrack(i3);
2416     phi3=fTrack->Phi();
2417     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2418     for(Int_t i4=0;i4<dMult;i4++)
2419     {
2420      if(i4==i1||i4==i2||i4==i3)continue;
2421      fTrack=anEvent->GetTrack(i4);
2422      phi4=fTrack->Phi();
2423      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2424      for(Int_t i5=0;i5<dMult;i5++)
2425      {
2426       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2427       fTrack=anEvent->GetTrack(i5);
2428       phi5=fTrack->Phi();
2429       if(phiWeights) wPhi5 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi5*nBinsPhi/TMath::TwoPi())));
2430       for(Int_t i6=0;i6<dMult;i6++)
2431       {
2432        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2433        fTrack=anEvent->GetTrack(i6);
2434        phi6=fTrack->Phi();
2435        if(phiWeights) wPhi6 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi6*nBinsPhi/TMath::TwoPi()))); 
2436        for(Int_t i7=0;i7<dMult;i7++)
2437        {
2438         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2439         fTrack=anEvent->GetTrack(i7);
2440         phi7=fTrack->Phi();
2441         if(phiWeights) wPhi7 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi7*nBinsPhi/TMath::TwoPi()))); 
2442         for(Int_t i8=0;i8<dMult;i8++)
2443         {
2444          if(i8==i1||i8==i2||i8==i3||i8==i4||i8==i5||i8==i6||i8==i7)continue;
2445          fTrack=anEvent->GetTrack(i8);
2446          phi8=fTrack->Phi();
2447          if(phiWeights) wPhi8 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi8*nBinsPhi/TMath::TwoPi())));  
2448          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}
2449         } 
2450        } 
2451       } 
2452      }
2453     }  
2454    }
2455   }
2456  }
2457  
2458  */
2459
2460  
2461
2462  // binning details of fDirectCorrelations (differential flow):
2463  //
2464  //101st bin: <2'>_{n|n}
2465  
2466  
2467  /*
2468  
2469  //<2'>_{n|n}
2470  for(Int_t i1=0;i1<nPrim;i1++)
2471  {
2472   fTrack=anEvent->GetTrack(i1);
2473   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2474   phi1=fTrack->Phi();
2475   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2476   for(Int_t i2=0;i2<nPrim;i2++)
2477   {
2478    if(i2==i1)continue;
2479    fTrack=anEvent->GetTrack(i2);
2480    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition
2481    phi2=fTrack->Phi();
2482    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2483     //cout<<"1st = "<<i1<<"     "<< (anEvent->GetTrack(i1))->Eta() << " " << (anEvent->GetTrack(i1))->Pt()<<endl;
2484     //cout<<"2nd = "<<i2<<"     "<< (anEvent->GetTrack(i2))->Eta() << " " << (anEvent->GetTrack(i2))->Pt()<<endl; 
2485    //fill the fDirectCorrelations:    
2486    fDirectCorrelations->Fill(100.,cos(1.*n*(phi1-phi2)),wPhi2);//<2'>_{n,n}
2487    fDirectCorrelations->Fill(103.,cos(1.*n*(phi1-phi2)),pow(wPhi1,2)*wPhi2);//<2'>_{n,n}
2488    fDirectCorrelations->Fill(104.,cos(2.*n*(phi1-phi2)),wPhi1*pow(wPhi2,2));//<2'>_{n,n}
2489    fDirectCorrelations->Fill(105.,cos(1.*n*(phi1-phi2)),pow(wPhi2,3));//<2'>_{n,n}
2490    
2491    
2492    
2493    fDirectCorrelations->Fill(41.,cos(2.*n*(phi1-phi2)),1);//<2'>_{2n,2n}
2494    fDirectCorrelations->Fill(42.,cos(3.*n*(phi1-phi2)),1);//<2'>_{3n,3n}
2495    fDirectCorrelations->Fill(43.,cos(4.*n*(phi1-phi2)),1);//<2'>_{4n,4n}   
2496     
2497   }//end of for(Int_t i2=0;i2<nPrim;i2++)
2498  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2499
2500  */ 
2501   
2502  
2503  
2504  /*
2505  
2506  //<3'>_{2n|n,n}
2507  for(Int_t i1=0;i1<nPrim;i1++)
2508  {
2509   fTrack=anEvent->GetTrack(i1);
2510   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2511   phi1=fTrack->Phi();
2512   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2513   for(Int_t i2=0;i2<nPrim;i2++)
2514   {
2515    if(i2==i1)continue;
2516    fTrack=anEvent->GetTrack(i2);
2517    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition
2518    phi2=fTrack->Phi();
2519    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2520    for(Int_t i3=0;i3<nPrim;i3++)
2521    {
2522     if(i3==i1||i3==i2)continue;
2523     fTrack=anEvent->GetTrack(i3);
2524     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition
2525     phi3=fTrack->Phi();
2526     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2527     //fill the fDirectCorrelations:     
2528     
2529     // 2-p
2530     fDirectCorrelations->Fill(101.,cos(n*(phi2-phi3)),wPhi1*wPhi2*wPhi3); // <w1 w2 w3 cos(n(phi2-phi3))>
2531     fDirectCorrelations->Fill(102.,cos(n*(phi1-phi3)),pow(wPhi2,2.)*wPhi3); // <w2^2 w3 cos(n(psi1-phi2))>
2532     
2533     // 3-p            
2534     fDirectCorrelations->Fill(110.,cos(n*(2.*phi1-phi2-phi3)),wPhi1*wPhi2*wPhi3); // <w1 w2 w3 cos(n(2psi1-phi2-phi3))>
2535     fDirectCorrelations->Fill(111.,cos(n*(phi1+phi2-2.*phi3)),wPhi2*pow(wPhi3,2.)); // <w2 w3^2 cos(n(psi1+phi2-2.*phi3))>
2536     
2537     
2538     //fDirectCorrelations->Fill(46.,cos(n*(phi1+phi2-2.*phi3)),1);//<3'>_{n,n|2n}    
2539    }//end of for(Int_t i3=0;i3<nPrim;i3++)  
2540   }//end of for(Int_t i2=0;i2<nPrim;i2++)  
2541  }//end of for(Int_t i1=0;i1<nPrim;i1++) 
2542  
2543  
2544  */
2545  
2546  
2547  /*
2548  
2549  //<4'>_{n,n|n,n}
2550  for(Int_t i1=0;i1<nPrim;i1++)
2551  {
2552   fTrack=anEvent->GetTrack(i1);
2553   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2554   tempCounter++;
2555   phi1=fTrack->Phi();
2556   if(phiWeights) wPhi1 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi1*nBinsPhi/TMath::TwoPi())));
2557   for(Int_t i2=0;i2<nPrim;i2++)
2558   {
2559    if(i2==i1)continue;
2560    fTrack=anEvent->GetTrack(i2);
2561    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2562    phi2=fTrack->Phi();
2563    if(phiWeights) wPhi2 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*nBinsPhi/TMath::TwoPi())));
2564    for(Int_t i3=0;i3<nPrim;i3++)
2565    { 
2566     if(i3==i1||i3==i2)continue;
2567     fTrack=anEvent->GetTrack(i3);
2568     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2569     phi3=fTrack->Phi();
2570     if(phiWeights) wPhi3 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*nBinsPhi/TMath::TwoPi())));
2571     for(Int_t i4=0;i4<nPrim;i4++)
2572     {
2573      if(i4==i1||i4==i2||i4==i3)continue;
2574      fTrack=anEvent->GetTrack(i4);
2575      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2576      phi4=fTrack->Phi();
2577      if(phiWeights) wPhi4 = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*nBinsPhi/TMath::TwoPi())));
2578      //fill the fDirectCorrelations:
2579      // 4-p            
2580      fDirectCorrelations->Fill(120.,cos(n*(phi1+phi2-phi3-phi4)),wPhi2*wPhi3*wPhi4); // <w1 w2 w3 cos(n(psi1+phi1-phi2-phi3))> 
2581      
2582        tempLoop+=wPhi2*wPhi3*wPhi4;
2583      
2584     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2585    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2586   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2587  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2588   
2589  */
2590  
2591  
2592  
2593  
2594  
2595  
2596  
2597  
2598  
2599  
2600  
2601  
2602  
2603  
2604  
2605  
2606  /*
2607  
2608  
2609  
2610  
2611  
2612  
2613  
2614  //<5'>_{2n,n|n,n,n}
2615  for(Int_t i1=0;i1<nPrim;i1++)
2616  {
2617   fTrack=anEvent->GetTrack(i1);
2618   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2619   phi1=fTrack->Phi();
2620   for(Int_t i2=0;i2<nPrim;i2++)
2621   {
2622    if(i2==i1)continue;
2623    fTrack=anEvent->GetTrack(i2);
2624    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2625    phi2=fTrack->Phi();
2626    for(Int_t i3=0;i3<nPrim;i3++)
2627    { 
2628     if(i3==i1||i3==i2)continue;
2629     fTrack=anEvent->GetTrack(i3);
2630     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2631     phi3=fTrack->Phi();
2632     for(Int_t i4=0;i4<nPrim;i4++)
2633     {
2634      if(i4==i1||i4==i2||i4==i3)continue;
2635      fTrack=anEvent->GetTrack(i4);
2636      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2637      phi4=fTrack->Phi();//
2638      for(Int_t i5=0;i5<nPrim;i5++)
2639      {
2640       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2641       fTrack=anEvent->GetTrack(i5);
2642       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2643       phi5=fTrack->Phi();    
2644       //fill the fDirectCorrelations:if(bNestedLoops)
2645       fDirectCorrelations->Fill(55.,cos(2.*n*phi1+n*phi2-n*phi3-n*phi4-n*phi5),1);//<5'>_{2n,n|n,n,n}
2646      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2647     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2648    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2649   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2650  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2651  
2652  //<6'>_{n,n,n|n,n,n}
2653  for(Int_t i1=0;i1<nPrim;i1++)
2654  {
2655   fTrack=anEvent->GetTrack(i1);
2656   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2657   phi1=fTrack->Phi();
2658   for(Int_t i2=0;i2<nPrim;i2++)
2659   {
2660    if(i2==i1)continue;
2661    fTrack=anEvent->GetTrack(i2);
2662    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2663    phi2=fTrack->Phi();
2664    for(Int_t i3=0;i3<nPrim;i3++)
2665    { 
2666     if(i3==i1||i3==i2)continue;
2667     fTrack=anEvent->GetTrack(i3);
2668     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2669     phi3=fTrack->Phi();
2670     for(Int_t i4=0;i4<nPrim;i4++)
2671     {
2672      if(i4==i1||i4==i2||i4==i3)continue;
2673      fTrack=anEvent->GetTrack(i4);
2674      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2675      phi4=fTrack->Phi();
2676      for(Int_t i5=0;i5<nPrim;i5++)
2677      {
2678       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2679       fTrack=anEvent->GetTrack(i5);
2680       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2681       phi5=fTrack->Phi();    
2682       for(Int_t i6=0;i6<nPrim;i6++)
2683       {
2684        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2685        fTrack=anEvent->GetTrack(i6);
2686        if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2687        phi6=fTrack->Phi();  
2688        //fill the fDirectCorrelations:
2689        fDirectCorrelations->Fill(60.,cos(n*(phi1+phi2+phi3-phi4-phi5-phi6)),1);//<6'>_{n,n,n|n,n,n}
2690       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
2691      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2692     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2693    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2694   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2695  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2696  
2697  //<7'>_{2n,n,n|n,n,n,n}
2698  for(Int_t i1=0;i1<nPrim;i1++)
2699  {
2700   fTrack=anEvent->GetTrack(i1);
2701   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2702   phi1=fTrack->Phi();
2703   for(Int_t i2=0;i2<nPrim;i2++)
2704   {
2705    if(i2==i1)continue;
2706    fTrack=anEvent->GetTrack(i2);
2707    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2708    phi2=fTrack->Phi();
2709    for(Int_t i3=0;i3<nPrim;i3++)
2710    { 
2711     if(i3==i1||i3==i2)continue;
2712     fTrack=anEvent->GetTrack(i3);
2713     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2714     phi3=fTrack->Phi();
2715     for(Int_t i4=0;i4<nPrim;i4++)
2716     {
2717      if(i4==i1||i4==i2||i4==i3)continue;
2718      fTrack=anEvent->GetTrack(i4);
2719      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2720      phi4=fTrack->Phi();
2721      for(Int_t i5=0;i5<nPrim;i5++)
2722      {
2723       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2724       fTrack=anEvent->GetTrack(i5);
2725       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2726       phi5=fTrack->Phi();    
2727       for(Int_t i6=0;i6<nPrim;i6++)
2728       {
2729        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2730        fTrack=anEvent->GetTrack(i6);
2731        if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2732        phi6=fTrack->Phi();
2733        for(Int_t i7=0;i7<nPrim;i7++)
2734        {
2735         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2736         fTrack=anEvent->GetTrack(i7);
2737         if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2738         phi7=fTrack->Phi();   
2739         //fill the fDirectCorrelations:
2740         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}
2741        }//end of for(Int_t i7=0;i7<nPrim;i7++)  
2742       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
2743      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2744     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2745    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2746   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2747  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2748  
2749  //<8'>_{n,n,n,n|n,n,n,n}
2750  for(Int_t i1=0;i1<nPrim;i1++)
2751  {
2752   fTrack=anEvent->GetTrack(i1);
2753   if(!((fTrack->Pt()>=0.5&&fTrack->Pt()<0.6)&&(fTrack->UseForDifferentialFlow())))continue;//POI condition
2754   phi1=fTrack->Phi();
2755   for(Int_t i2=0;i2<nPrim;i2++)
2756   {
2757    if(i2==i1)continue;
2758    fTrack=anEvent->GetTrack(i2);
2759    if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2760    phi2=fTrack->Phi();
2761    for(Int_t i3=0;i3<nPrim;i3++)
2762    { 
2763     if(i3==i1||i3==i2)continue;
2764     fTrack=anEvent->GetTrack(i3);
2765     if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition   
2766     phi3=fTrack->Phi();
2767     for(Int_t i4=0;i4<nPrim;i4++)
2768     {
2769      if(i4==i1||i4==i2||i4==i3)continue;
2770      fTrack=anEvent->GetTrack(i4);
2771      if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2772      phi4=fTrack->Phi();
2773      for(Int_t i5=0;i5<nPrim;i5++)
2774      {
2775       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
2776       fTrack=anEvent->GetTrack(i5);
2777       if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2778       phi5=fTrack->Phi();    
2779       for(Int_t i6=0;i6<nPrim;i6++)
2780       {
2781        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
2782        fTrack=anEvent->GetTrack(i6);
2783        if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2784        phi6=fTrack->Phi();
2785        for(Int_t i7=0;i7<nPrim;i7++)
2786        {
2787         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
2788         fTrack=anEvent->GetTrack(i7);
2789         if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2790         phi7=fTrack->Phi();
2791         for(Int_t i8=0;i8<nPrim;i8++)
2792         {
2793          if(i8==i1||i8==i2||i8==i3||i8==i4||i8==i5||i8==i6||i8==i7)continue;
2794          fTrack=anEvent->GetTrack(i8);
2795          if(!(fTrack->UseForIntegratedFlow()))continue;//RP condition  
2796          phi8=fTrack->Phi();           
2797          //fill the fDirectCorrelations:
2798          fDirectCorrelations->Fill(70.,cos(n*(phi1+phi2+phi3+phi4-phi5-phi6-phi7-phi8)),1);//<8'>_{n,n,n,n|n,n,n,n}
2799         }//end of for(Int_t i8=0;i8<nPrim;i8++) 
2800        }//end of for(Int_t i7=0;i7<nPrim;i7++)  
2801       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
2802      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
2803     }//end of for(Int_t i4=0;i4<nPrim;i4++)
2804    }//end of for(Int_t i3=0;i3<nPrim;i3++)
2805   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
2806  }//end of for(Int_t i1=0;i1<nPrim;i1++)
2807  
2808  
2809  
2810  
2811  
2812  
2813  
2814  
2815  
2816  
2817  */
2818  
2819  
2820  
2821  
2822  
2823  
2824  
2825  
2826  
2827  
2828  
2829  
2830  
2831  
2832  
2833  
2834  
2835  
2836  
2837  
2838  
2839  
2840
2841  
2842  
2843  
2844  
2845  
2846  
2847  
2848  
2849  
2850 //--------------------------------------------------------------------------------------------------------------------------------  
2851
2852  //}//end of if(nPrim>0&&nPrim<12)
2853 }//end of Make()
2854
2855 //================================================================================================================
2856
2857 void AliFlowAnalysisWithQCumulants::Finish()
2858 {
2859  //calculate the final results
2860  
2861  // harmonics
2862  Int_t n = 2; // to be improved 
2863  
2864  //--------------------------------------------------------------------------------------------------------- 
2865  // avarage multiplicity
2866  Double_t AvMPOI = (fCommonHists2nd->GetHistMultDiff())->GetMean(); // to be improved 
2867  Double_t AvMRP  = (fCommonHists2nd->GetHistMultInt())->GetMean(); // to be improved 
2868
2869  // number of events
2870  Double_t nEvtsPOI = (fCommonHists2nd->GetHistMultDiff())->GetEntries(); // to be improved
2871  Double_t nEvtsRP  = (fCommonHists2nd->GetHistMultInt())->GetEntries(); // to be improved
2872  //---------------------------------------------------------------------------------------------------------
2873  
2874  //---------------------------------------------------------------------------------------------------------
2875  // 2-, 4-, 6- and 8-particle azimuthal correlation:
2876  Double_t two   = fQCorrelations->GetBinContent(1);  //<<2>>_{n|n}
2877  Double_t four  = fQCorrelations->GetBinContent(11); //<<4>>_{n,n|n,n}
2878  Double_t six   = fQCorrelations->GetBinContent(24); //<<6>>_{n,n,n|n,n,n}
2879  Double_t eight = fQCorrelations->GetBinContent(31); //<<8>>_{n,n,n,n|n,n,n,n}
2880  
2881  // 2nd, 4th, 6th and 8th order Q-cumulant:
2882  Double_t secondOrderQCumulant = two; //c_n{2} 
2883  Double_t fourthOrderQCumulant = four-2.*pow(two,2.); //c_n{4}
2884  Double_t sixthOrderQCumulant  = six-9.*two*four+12.*pow(two,3.); //c_n{6}
2885  Double_t eightOrderQCumulant  = eight-16.*two*six-18.*pow(four,2.)+144.*pow(two,2.)*four-144.*pow(two,4.); //c_n{8} 
2886  
2887  // "no-name" integrated flow estimates from Q-cumulants:
2888  Double_t vn2=0.,vn4=0.,vn6=0.,vn8=0.;
2889  // Double_t sd2=0.,sd4=0.,sd6=0.,sd8=0.; 
2890  Double_t sd6=0.,sd8=0.; // to be improved/removed
2891  if(secondOrderQCumulant>0.)
2892  {
2893   vn2 = pow(secondOrderQCumulant,0.5); //v_n{2}
2894  } 
2895  if(fourthOrderQCumulant<0.)
2896  {
2897   vn4 = pow(-fourthOrderQCumulant,1./4.); //v_n{4}
2898  } 
2899  if(sixthOrderQCumulant>0.)
2900  {
2901   vn6 = pow((1./4.)*sixthOrderQCumulant,1./6.); //v_n{6}
2902  } 
2903  if(eightOrderQCumulant<0.)
2904  {
2905   vn8 = pow((-1./33.)*eightOrderQCumulant,1./8.); //v_n{8}
2906  }
2907  //---------------------------------------------------------------------------------------------------------
2908
2909  //--------------------------------------------------------------------------------------------------------- 
2910  // weighted 2-, 4-, 6- and 8-particle azimuthal correlation:
2911  Double_t twoW   = fWeightedQCorrelations->GetBinContent(1);  //<<2>>_{n|n}
2912  Double_t fourW  = fWeightedQCorrelations->GetBinContent(21); //<<4>>_{n,n|n,n}
2913  // Double_t sixW   = fWeightedQCorrelations->GetBinContent(24); //<<6>>_{n,n,n|n,n,n}
2914  // Double_t eightW = fWeightedQCorrelations->GetBinContent(31); //<<8>>_{n,n,n,n|n,n,n,n}
2915  
2916  // 2nd, 4th, 6th and 8th order weighted Q-cumulant:
2917  Double_t secondOrderQCumulantW = twoW; //c_n{2} 
2918  Double_t fourthOrderQCumulantW = fourW-2.*pow(twoW,2.); //c_n{4}
2919  // Double_t sixthOrderQCumulantW  = sixW-9.*twoW*fourW+12.*pow(twoW,3.); //c_n{6}
2920  // Double_t eightOrderQCumulantW  = eightW-16.*twoW*sixW-18.*pow(fourW,2.)+144.*pow(twoW,2.)*fourW-144.*pow(twoW,4.); //c_n{8} 
2921  
2922  // "no-name" integrated flow estimates from weighted Q-cumulants:
2923  cout<<endl;
2924  cout<<"**************************************"<<endl;
2925  cout<<"**************************************"<<endl;
2926  cout<<"flow estimates from Q-cumulants :"<<endl;
2927  cout<<endl;
2928  
2929  Double_t vn2W=0.,vn4W=0.;
2930  Double_t sd2W=0.,sd4W=0.; 
2931  if(secondOrderQCumulantW>0.){
2932   vn2W = pow(secondOrderQCumulantW,0.5); // weighted v_n{2}
2933   //sd2W = 0.5*pow(secondOrderQCumulantW,-0.5)*secondOrderQCumulantErrorW; // to be improved (correct treatment of errors needed)
2934   cout<<" v_"<<n<<"{2} = "<<vn2W<<" +/- "<<sd2W<<endl;
2935   fIntFlowResultsQC->SetBinContent(1,vn2W);
2936   //fIntFlowResultsQC->SetBinError(1,sd2W);
2937   //common histograms:
2938   fCommonHistsResults2nd->FillIntegratedFlow(vn2W,sd2W);
2939   //fCommonHistsResults2nd->FillChi(vn2W*pow(AvM,0.5)); // to be removed
2940  }else{
2941   cout<<" v_"<<n<<"{2} = Im"<<endl;
2942  }          
2943  if(fourW!=0. && fourthOrderQCumulantW<0.){
2944   vn4W = pow(-fourthOrderQCumulantW,1./4.); //v_n{4}
2945   //sd4W = 0.25*pow(-fourthOrderQCumulantW,-3./4.)*fourthOrderQCumulantErrorW; // to be improved (correct treatment of errors needed)
2946   cout<<" v_"<<n<<"{4} = "<<vn4W<<" +/- "<<sd4W<<endl;
2947   fIntFlowResultsQC->SetBinContent(2,vn4W);
2948   //fIntFlowResultsQC->SetBinError(2,sd4W);
2949   //common histograms:
2950   fCommonHistsResults4th->FillIntegratedFlow(vn4W,sd4W);
2951   //fCommonHistsResults4th->FillChi(vn4W*pow(AvM,0.5)); // to be removed
2952  }else{
2953   cout<<" v_"<<n<<"{4} = Im"<<endl;
2954  }
2955  // !!! to be improved (6th and 8th order are without weights) !!!
2956  if(six!=0. && sixthOrderQCumulant>0.){
2957   vn6 = pow((1./4.)*sixthOrderQCumulant,1./6.); //v_n{6}
2958   //sd6 = (1./6.)*pow(2.,-1./3.)*pow(sixthOrderQCumulant,-5./6.)*sixthOrderQCumulantError;
2959   cout<<" v_"<<n<<"{6} = "<<vn6<<" +/- "<<sd6<<endl;
2960   fIntFlowResultsQC->SetBinContent(3,vn6);
2961   //fIntFlowResultsQC->SetBinError(3,sd6);
2962   //common histograms:
2963   fCommonHistsResults6th->FillIntegratedFlow(vn6,sd6);
2964   //fCommonHistsResults6th->FillChi(vn6*pow(AvM,0.5));//to be removed
2965  }else{
2966   cout<<" v_"<<n<<"{6} = Im"<<endl;
2967  }
2968  if(eight!=0. && eightOrderQCumulant<0.){
2969   vn8 = pow((-1./33.)*eightOrderQCumulant,1./8.); //v_n{8}
2970   cout<<" v_"<<n<<"{8} = "<<vn8<<" +/- "<<sd8<<endl;
2971   fIntFlowResultsQC->SetBinContent(4,vn8);
2972   //fIntFlowResultsQC->SetBinError(4,sd8);
2973   //common histograms:
2974   fCommonHistsResults8th->FillIntegratedFlow(vn8,sd8);
2975   //fCommonHistsResults8th->FillChi(vn8*pow(AvM,0.5));//to be removed
2976  }else{
2977   cout<<" v_"<<n<<"{8} = Im"<<endl;
2978  }
2979  cout<<endl;
2980  cout<<"   nEvts = "<<nEvtsRP<<", AvM = "<<AvMRP<<endl; // to be improved
2981  cout<<"**************************************"<<endl;
2982  cout<<"**************************************"<<endl;
2983  cout<<endl; 
2984 //--------------------------------------------------------------------------------------------------------- 
2985  
2986 //---------------------------------------------------------------------------------------------------------
2987 // differential flow (POI)
2988 Int_t nBinsPtPOI = f2WPerPtBin1n1nPOI->GetNbinsX();
2989 Int_t nBinsEtaPOI = f4WPerPtBin1n1n1n1nPOI->GetNbinsX();
2990
2991 // Pt:
2992 Double_t secondOrderQCumulantDiffFlowPtPOIW = 0.;
2993 Double_t fourthOrderQCumulantDiffFlowPtPOIW = 0.;
2994
2995 Double_t dVn2ndPOIW=0.,dSd2ndPOIW=0.,dDiffvn2ndPOIW=0.,dYield2ndPOIW=0.,dSum2ndPOIW=0.;
2996 Double_t dVn4thPOIW=0.,dSd4thPOIW=0.,dDiffvn4thPOIW=0.,dYield4thPOIW=0.,dSum4thPOIW=0.;
2997
2998 for(Int_t bb=1;bb<nBinsPtPOI+1;bb++)
2999 {
3000  // QC{2}
3001  if(f2WPerPtBin1n1nPOI->GetBinEntries(bb)>0.&&vn2W!=0)
3002  { 
3003   secondOrderQCumulantDiffFlowPtPOIW = f2WPerPtBin1n1nPOI->GetBinContent(bb); // with weights
3004   fDiffFlowResults2ndOrderQC->SetBinContent(bb,secondOrderQCumulantDiffFlowPtPOIW/vn2W);
3005   // common histogram:
3006   fCommonHistsResults2nd->FillDifferentialFlowPtPOI(bb,secondOrderQCumulantDiffFlowPtPOIW/vn2W, 0.); //to be improved (errors && bb or bb+1 ?)
3007   // -------------------------------------------------------------------
3008   // integrated flow (weighted, POI, Pt, 2nd order):
3009   dDiffvn2ndPOIW=(fCommonHistsResults2nd->GetHistDiffFlowPtPOI())->GetBinContent(bb);
3010   dYield2ndPOIW=(fCommonHists2nd->GetHistPtDiff())->GetBinContent(bb);
3011   dVn2ndPOIW+=dDiffvn2ndPOIW*dYield2ndPOIW;
3012   dSum2ndPOIW+=dYield2ndPOIW;
3013   // -------------------------------------------------------------------
3014  }
3015  // QC{4]
3016  if(f4WPerPtBin1n1n1n1nPOI->GetBinEntries(bb)>0.&&vn4W!=0.)
3017  {
3018   fourthOrderQCumulantDiffFlowPtPOIW = f4WPerPtBin1n1n1n1nPOI->GetBinContent(bb)-2.*f2WPerPtBin1n1nPOI->GetBinContent(bb)*pow(vn2W,2.); // with weights
3019   fDiffFlowResults4thOrderQC->SetBinContent(bb,-1.*fourthOrderQCumulantDiffFlowPtPOIW/pow(vn4W,3.));
3020   //common histogram:
3021   fCommonHistsResults4th->FillDifferentialFlowPtPOI(bb,-1.*fourthOrderQCumulantDiffFlowPtPOIW/pow(vn4W,3.), 0.); //to be improved (errors)
3022   // -------------------------------------------------------------------
3023   //integrated flow (POI, Pt, 4th order):
3024   dDiffvn4thPOIW=(fCommonHistsResults4th->GetHistDiffFlowPtPOI())->GetBinContent(bb);
3025   dYield4thPOIW=(fCommonHists4th->GetHistPtDiff())->GetBinContent(bb);
3026   dVn4thPOIW+=dDiffvn4thPOIW*dYield4thPOIW;
3027   dSum4thPOIW+=dYield4thPOIW;
3028   // -------------------------------------------------------------------
3029  }
3030 }      
3031
3032 cout<<endl;
3033 cout<<"**************************************"<<endl;
3034 cout<<"**************************************"<<endl;
3035 cout<<"flow estimates from Q-cumulants (POI):"<<endl;
3036 cout<<endl;
3037 //storing the final results for integrated flow (POI):
3038 // QC{2}
3039 if(dSum2ndPOIW && fCommonHistsResults2nd)
3040 {
3041  dVn2ndPOIW/=dSum2ndPOIW;
3042  fCommonHistsResults2nd->FillIntegratedFlowPOI(dVn2ndPOIW,0.); // to be improved (errors)
3043  cout<<" v_"<<n<<"{2} = "<<dVn2ndPOIW<<" +/- "<<dSd2ndPOIW<<endl;
3044 }else 
3045  {
3046   cout<<" v_"<<n<<"{2} = Im"<<endl;
3047  }
3048
3049 // QC{4}
3050 if(dSum4thPOIW && fCommonHistsResults4th)
3051 {
3052  dVn4thPOIW/=dSum4thPOIW;
3053  fCommonHistsResults4th->FillIntegratedFlowPOI(dVn4thPOIW,0.); // to be improved (errors)
3054  cout<<" v_"<<n<<"{4} = "<<dVn4thPOIW<<" +/- "<<dSd4thPOIW<<endl;
3055 }else
3056  {
3057   cout<<" v_"<<n<<"{4} = Im"<<endl;
3058  }
3059
3060 cout<<endl;
3061 cout<<"   nEvts = "<<nEvtsPOI<<", AvM = "<<AvMPOI<<endl;
3062 cout<<"**************************************"<<endl;
3063 cout<<"**************************************"<<endl;
3064 cout<<endl;  
3065
3066 //Eta:
3067 Double_t secondOrderQCumulantDiffFlowEtaPOIW = 0.;
3068 Double_t fourthOrderQCumulantDiffFlowEtaPOIW = 0.;
3069
3070 for(Int_t bb=1;bb<nBinsEtaPOI+1;bb++)
3071 {
3072  if(f2WPerEtaBin1n1nPOI->GetBinEntries(bb)>0.&&vn2W!=0)
3073  {
3074   secondOrderQCumulantDiffFlowEtaPOIW = f2WPerEtaBin1n1nPOI->GetBinContent(bb); // with weights
3075   fDiffFlowResults2ndOrderQC->SetBinContent(bb,secondOrderQCumulantDiffFlowEtaPOIW/vn2W);
3076   //common histogram:
3077   fCommonHistsResults2nd->FillDifferentialFlowEtaPOI(bb,secondOrderQCumulantDiffFlowEtaPOIW/vn2W, 0.);//to be improved (errors)
3078  }
3079  if(f4WPerEtaBin1n1n1n1nPOI->GetBinEntries(bb)>0.&&vn4W!=0.)
3080  {
3081   fourthOrderQCumulantDiffFlowEtaPOIW = f4WPerEtaBin1n1n1n1nPOI->GetBinContent(bb)-2.*f2WPerEtaBin1n1nPOI->GetBinContent(bb)*pow(vn2W,2.); // with weights
3082   fDiffFlowResults4thOrderQC->SetBinContent(bb,-1.*fourthOrderQCumulantDiffFlowEtaPOIW/pow(vn4W,3.));
3083   //common histogram:
3084   fCommonHistsResults4th->FillDifferentialFlowEtaPOI(bb,-1.*fourthOrderQCumulantDiffFlowEtaPOIW/pow(vn4W,3.), 0.);//to be improved (errors)
3085  }
3086 }      
3087 //------------------------------------------------------------
3088  
3089  
3090  
3091  
3092  
3093  //---------------------------------------------------------------------------------------------------------
3094 // differential flow (RP)
3095 Int_t nBinsPtRP = f2WPerPtBin1n1nRP->GetNbinsX();
3096 Int_t nBinsEtaRP = f4WPerPtBin1n1n1n1nRP->GetNbinsX();
3097
3098 // Pt:
3099 Double_t secondOrderQCumulantDiffFlowPtRPW = 0.;
3100 Double_t fourthOrderQCumulantDiffFlowPtRPW = 0.;
3101
3102 Double_t dVn2ndRPW=0.,dSd2ndRPW=0.,dDiffvn2ndRPW=0.,dYield2ndRPW=0.,dSum2ndRPW=0.;
3103 Double_t dVn4thRPW=0.,dSd4thRPW=0.,dDiffvn4thRPW=0.,dYield4thRPW=0.,dSum4thRPW=0.;
3104
3105 for(Int_t bb=1;bb<nBinsPtRP+1;bb++)
3106 {
3107  // QC{2}
3108  if(f2WPerPtBin1n1nRP->GetBinEntries(bb)>0.&&vn2W!=0)
3109  { 
3110   secondOrderQCumulantDiffFlowPtRPW = f2WPerPtBin1n1nRP->GetBinContent(bb); // with weights
3111   fDiffFlowResults2ndOrderQC->SetBinContent(bb,secondOrderQCumulantDiffFlowPtRPW/vn2W);
3112   // common histogram:
3113   fCommonHistsResults2nd->FillDifferentialFlowPtRP(bb,secondOrderQCumulantDiffFlowPtRPW/vn2W, 0.); //to be improved (errors && bb or bb+1 ?)
3114   // -------------------------------------------------------------------
3115   // integrated flow (weighted, RP, Pt, 2nd order):
3116   dDiffvn2ndRPW=(fCommonHistsResults2nd->GetHistDiffFlowPtRP())->GetBinContent(bb);
3117   dYield2ndRPW=(fCommonHists2nd->GetHistPtInt())->GetBinContent(bb);
3118   dVn2ndRPW+=dDiffvn2ndRPW*dYield2ndRPW;
3119   dSum2ndRPW+=dYield2ndRPW;
3120   // -------------------------------------------------------------------
3121  }
3122  // QC{4]
3123  if(f4WPerPtBin1n1n1n1nRP->GetBinEntries(bb)>0.&&vn4W!=0.)
3124  {
3125   fourthOrderQCumulantDiffFlowPtRPW = f4WPerPtBin1n1n1n1nRP->GetBinContent(bb)-2.*f2WPerPtBin1n1nRP->GetBinContent(bb)*pow(vn2W,2.); // with weights
3126   fDiffFlowResults4thOrderQC->SetBinContent(bb,-1.*fourthOrderQCumulantDiffFlowPtRPW/pow(vn4W,3.));
3127   //common histogram:
3128   fCommonHistsResults4th->FillDifferentialFlowPtRP(bb,-1.*fourthOrderQCumulantDiffFlowPtRPW/pow(vn4W,3.), 0.); //to be improved (errors)
3129   // -------------------------------------------------------------------
3130   //integrated flow (RP, Pt, 4th order):
3131   dDiffvn4thRPW=(fCommonHistsResults4th->GetHistDiffFlowPtRP())->GetBinContent(bb);
3132   dYield4thRPW=(fCommonHists4th->GetHistPtInt())->GetBinContent(bb);
3133   dVn4thRPW+=dDiffvn4thRPW*dYield4thRPW;
3134   dSum4thRPW+=dYield4thRPW;
3135   // -------------------------------------------------------------------
3136  }
3137 }      
3138
3139 cout<<endl;
3140 cout<<"**************************************"<<endl;
3141 cout<<"**************************************"<<endl;
3142 cout<<"flow estimates from Q-cumulants (RP):"<<endl;
3143 cout<<endl;
3144 //storing the final results for integrated flow (RP):
3145 // QC{2}
3146 if(dSum2ndRPW && fCommonHistsResults2nd)
3147 {
3148  dVn2ndRPW/=dSum2ndRPW;
3149  fCommonHistsResults2nd->FillIntegratedFlowRP(dVn2ndRPW,0.); // to be improved (errors)
3150  cout<<" v_"<<n<<"{2} = "<<dVn2ndRPW<<" +/- "<<dSd2ndRPW<<endl;
3151 }else 
3152  {
3153   cout<<" v_"<<n<<"{2} = Im"<<endl;
3154  }