fix formatting for microsoft compiler
[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
31 #include "TFile.h"
32 #include "TList.h"
33 #include "TGraph.h"
34 #include "TParticle.h"
35 #include "TRandom3.h"
36 #include "TStyle.h"
37 #include "TProfile.h"
38 #include "TProfile2D.h" 
39 #include "TProfile3D.h"
40 #include "TMath.h"
41 #include "TArrow.h"
42 #include "TPaveLabel.h"
43 #include "TCanvas.h"
44 #include "AliFlowEventSimple.h"
45 #include "AliFlowTrackSimple.h"
46 #include "AliFlowAnalysisWithQCumulants.h"
47 #include "TArrayD.h"
48 #include "TRandom.h"
49 #include "TF1.h"
50
51 class TH1;
52 class TH2;
53 class TGraph;
54 class TPave;
55 class TLatex;
56 class TMarker;
57 class TRandom3;
58 class TObjArray;
59 class TList;
60 class TCanvas;
61 class TSystem;
62 class TROOT;
63 class AliFlowVector;
64 class TVector;
65
66
67 //================================================================================================================
68
69
70 ClassImp(AliFlowAnalysisWithQCumulants)
71
72 AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants(): 
73  // 0.) base:
74  fHistList(NULL),
75  // 1.) common:
76  fCommonHists(NULL),
77  fCommonHists2nd(NULL), 
78  fCommonHists4th(NULL),
79  fCommonHists6th(NULL),
80  fCommonHists8th(NULL),
81  fCommonHistsResults2nd(NULL),
82  fCommonHistsResults4th(NULL),
83  fCommonHistsResults6th(NULL),
84  fCommonHistsResults8th(NULL),
85  fnBinsPhi(0),
86  fPhiMin(0),
87  fPhiMax(0),
88  fPhiBinWidth(0),
89  fnBinsPt(0),
90  fPtMin(0),
91  fPtMax(0),
92  fPtBinWidth(0),
93  fnBinsEta(0),
94  fEtaMin(0),
95  fEtaMax(0),
96  fEtaBinWidth(0),
97  fHarmonic(2),
98  fAnalysisLabel(NULL),
99  // 2.) weights:
100  fWeightsList(NULL),
101  fUsePhiWeights(kFALSE),
102  fUsePtWeights(kFALSE),
103  fUseEtaWeights(kFALSE),
104  fUseParticleWeights(NULL),
105  fPhiWeights(NULL),
106  fPtWeights(NULL),
107  fEtaWeights(NULL),
108  // 3.) integrated flow:
109  fIntFlowList(NULL), 
110  fIntFlowProfiles(NULL),
111  fIntFlowResults(NULL),
112  fIntFlowFlags(NULL),
113  fApplyCorrectionForNUA(kTRUE), 
114  fReQ(NULL),
115  fImQ(NULL),
116  fSMpk(NULL),
117  fIntFlowCorrelationsEBE(NULL),
118  fIntFlowEventWeightsForCorrelationsEBE(NULL),
119  fIntFlowCorrelationsAllEBE(NULL),
120  fAvMultiplicity(NULL),
121  fIntFlowCorrelationsPro(NULL),
122  fIntFlowCorrelationsAllPro(NULL),
123  fIntFlowExtraCorrelationsPro(NULL),
124  fIntFlowProductOfCorrelationsPro(NULL),
125  fIntFlowCorrelationsHist(NULL),
126  fIntFlowCorrelationsAllHist(NULL),
127  fIntFlowCovariances(NULL),
128  fIntFlowSumOfProductOfEventWeights(NULL),
129  fIntFlowQcumulants(NULL),
130  fIntFlow(NULL),
131  // 4.) differential flow:
132  fDiffFlowList(NULL),
133  fDiffFlowProfiles(NULL),
134  fDiffFlowResults(NULL),
135  fDiffFlowFlags(NULL),
136  fCalculate2DFlow(kFALSE),
137  // 5.) distributions:
138  fDistributionsList(NULL),
139  // x.) debugging and cross-checking:
140  fNestedLoopsList(NULL),
141  fEvaluateIntFlowNestedLoops(kFALSE),
142  fEvaluateDiffFlowNestedLoops(kFALSE),
143  fMaxAllowedMultiplicity(10),
144  fEvaluateNestedLoops(NULL),
145  fIntFlowDirectCorrelations(NULL),
146  fIntFlowExtraDirectCorrelations(NULL),
147  fCrossCheckInPtBinNo(10),
148  fCrossCheckInEtaBinNo(20)
149  {
150   // constructor  
151   
152   // base list to hold all output objects:
153   fHistList = new TList();
154   fHistList->SetName("cobjQC");
155   fHistList->SetOwner(kTRUE);
156   
157   // list to hold histograms with phi, pt and eta weights:      
158   fWeightsList = new TList();
159     
160   // analysis label;
161   fAnalysisLabel = new TString();
162       
163   // initialize all arrays:  
164   this->InitializeArraysForIntFlow();
165   this->InitializeArraysForDiffFlow();
166   this->InitializeArraysForDistributions();
167   this->InitializeArraysForNestedLoops();
168   
169  } // end of constructor
170  
171
172 //================================================================================================================  
173
174
175 AliFlowAnalysisWithQCumulants::~AliFlowAnalysisWithQCumulants()
176 {
177  // destructor
178  
179  delete fHistList;
180
181 } // end of AliFlowAnalysisWithQCumulants::~AliFlowAnalysisWithQCumulants()
182
183
184 //================================================================================================================
185
186
187 void AliFlowAnalysisWithQCumulants::Init()
188 {
189  // a) Access all common constants;
190  // b) Book all objects;
191  // c) Store flags for integrated and differential flow;
192  // d) Store harmonic which will be estimated.
193   
194  // a) Access all common constants:
195  this->AccessConstants();
196  
197  // b) Book all objects:
198  this->BookAndFillWeightsHistograms();
199  this->BookAndNestAllLists();
200  this->BookCommonHistograms();
201  this->BookEverythingForIntegratedFlow(); 
202  this->BookEverythingForDifferentialFlow(); 
203  this->BookEverythingForDistributions();
204  this->BookEverythingForNestedLoops();
205  
206  // c) Store flags for integrated and differential flow:
207  this->StoreIntFlowFlags();
208  this->StoreDiffFlowFlags();
209
210  // d) Store harmonic which will be estimated:
211  this->StoreHarmonic();
212  
213 } // end of void AliFlowAnalysisWithQCumulants::Init()
214
215
216 //================================================================================================================
217
218
219 void AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
220 {
221  // Running over data only in this method.
222  
223  //  a) Fill the common control histograms and call the method to fill fAvMultiplicity;
224  //  b) Loop over data and calculate e-b-e quantities;
225  //  c) call the methods;
226  //  d) Debugging and cross-checking (evaluate nested loops);
227  //  e) Reset all event by event quantities. 
228  
229  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
230  Double_t dPt  = 0.; // transverse momentum
231  Double_t dEta = 0.; // pseudorapidity
232
233  Double_t wPhi = 1.; // phi weight
234  Double_t wPt  = 1.; // pt weight
235  Double_t wEta = 1.; // eta weight
236  
237  Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of RPs (i.e. number of particles used to determine the reaction plane)
238  
239  // a) Fill the common control histograms and call the method to fill fAvMultiplicity:
240  this->FillCommonControlHistograms(anEvent);                                                               
241  this->FillAverageMultiplicities(nRP);                                                                  
242                                                                                                                                                                                                                                                                                         
243  // b) Loop over data and calculate e-b-e quantities:
244  Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where:
245                                            // nRP   = # of particles used to determine the reaction plane;
246                                            // nPOI  = # of particles of interest for a detailed flow analysis;
247                                            // rest  = # of particles which are not niether RPs nor POIs.  
248  
249  AliFlowTrackSimple *aftsTrack = NULL;
250  
251  for(Int_t i=0;i<nPrim;i++) 
252  { 
253   aftsTrack=anEvent->GetTrack(i);
254   if(aftsTrack)
255   {
256    if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())) continue; // consider only tracks which are RPs or POIs
257    Int_t n = fHarmonic; // shortcut for the harmonic
258    if(aftsTrack->InRPSelection()) // RP condition:
259    {    
260     dPhi = aftsTrack->Phi();
261     dPt  = aftsTrack->Pt();
262     dEta = aftsTrack->Eta();
263     if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi weight for this particle:
264     {
265      wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
266     }
267     if(fUsePtWeights && fPtWeights && fnBinsPt) // determine pt weight for this particle:
268     {
269      wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); 
270     }              
271     if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta weight for this particle: 
272     {
273      wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); 
274     } 
275       
276     // integrated flow: 
277     // calculate Re[Q_{m*n,k}] and Im[Q_{m*n,k}], m = 1,2,3,4, for this event:
278     for(Int_t m=0;m<4;m++)
279     {
280      for(Int_t k=0;k<9;k++)
281      {
282       (*fReQ)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1)*n*dPhi); 
283       (*fImQ)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1)*n*dPhi); 
284      } 
285     }
286     // calculate S^{M}_{p,k} for this event 
287     // Remark: final calculation of S^{M}_{p,k} follows after the loop over data bellow:
288     for(Int_t p=0;p<8;p++)
289     {
290      for(Int_t k=0;k<9;k++)
291      {     
292       (*fSMpk)(p,k)+=pow(wPhi*wPt*wEta,k);
293      }
294     } 
295     
296     // differential flow:
297     // 1D (pt):
298     // (r_{m*m,k}(pt)): 
299     for(Int_t m=0;m<4;m++)
300     {
301      for(Int_t k=0;k<9;k++)
302      {
303       fReRPQ1dEBE[0][0][m][k]->Fill(dPt,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
304       fImRPQ1dEBE[0][0][m][k]->Fill(dPt,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
305      }
306     }
307            
308     // s_{k}(pt) for RPs // to be improved (clarified)
309     // Remark: final calculation of s_{p,k}(pt) follows after the loop over data bellow:
310     for(Int_t k=0;k<9;k++)
311     {
312      fs1dEBE[0][0][k]->Fill(dPt,pow(wPhi*wPt*wEta,k),1.);
313     }
314     // 1D (eta):
315     // (r_{m*m,k}(eta)): 
316     for(Int_t m=0;m<4;m++)
317     {
318      for(Int_t k=0;k<9;k++)
319      {
320       fReRPQ1dEBE[0][1][m][k]->Fill(dEta,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
321       fImRPQ1dEBE[0][1][m][k]->Fill(dEta,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
322      }
323     }   
324     // s_{k}(eta) for RPs // to be improved (clarified)
325     // Remark: final calculation of s_{p,k}(eta) follows after the loop over data bellow:
326     for(Int_t k=0;k<9;k++)
327     {
328      fs1dEBE[0][1][k]->Fill(dEta,pow(wPhi*wPt*wEta,k),1.);
329     }
330     
331     
332     
333     /*
334     // 2D (pt,eta):
335     if(fCalculate2DFlow)
336     {
337      // (r_{m*m,k}(pt,eta)): 
338      for(Int_t m=0;m<4;m++)
339      {
340       for(Int_t k=0;k<9;k++)
341       {
342        fReRPQ2dEBE[0][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
343        fImRPQ2dEBE[0][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
344       }
345      }    
346      // s_{k}(pt,eta) for RPs // to be improved (clarified)
347      // Remark: final calculation of s_{p,k}(pt,eta) follows after the loop over data bellow:
348      for(Int_t k=0;k<9;k++)
349      {
350       fs2dEBE[0][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k),1.);
351      }
352     } // end of if(fCalculate2DFlow)  
353     */ 
354     
355       
356      
357     if(aftsTrack->InPOISelection())
358     {
359      // 1D (pt): 
360      // (q_{m*m,k}(pt)): 
361      for(Int_t m=0;m<4;m++)
362      {
363       for(Int_t k=0;k<9;k++)
364       {
365        fReRPQ1dEBE[2][0][m][k]->Fill(dPt,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
366        fImRPQ1dEBE[2][0][m][k]->Fill(dPt,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
367       }
368      } 
369      // s_{k}(pt) for RP&&POIs // to be improved (clarified)
370      // Remark: final calculation of s_{p,k}(pt,eta) follows after the loop over data bellow:
371      for(Int_t k=0;k<9;k++)
372      {
373       fs1dEBE[2][0][k]->Fill(dPt,pow(wPhi*wPt*wEta,k),1.);
374      }
375      // 1D (eta): 
376      // (q_{m*m,k}(eta)): 
377      for(Int_t m=0;m<4;m++)
378      {
379       for(Int_t k=0;k<9;k++)
380       {
381        fReRPQ1dEBE[2][1][m][k]->Fill(dEta,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
382        fImRPQ1dEBE[2][1][m][k]->Fill(dEta,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
383       }
384      } 
385      // s_{k}(eta) for RP&&POIs // to be improved (clarified)
386      // Remark: final calculation of s_{p,k}(pt,eta) follows after the loop over data bellow:
387      for(Int_t k=0;k<9;k++)
388      {
389       fs1dEBE[2][1][k]->Fill(dEta,pow(wPhi*wPt*wEta,k),1.);
390      }
391      
392      /*
393      // 2D (pt,eta) 
394      if(fCalculate2DFlow)
395      {
396       // (q_{m*m,k}(pt,eta)): 
397       for(Int_t m=0;m<4;m++)
398       {
399        for(Int_t k=0;k<9;k++)
400        {
401         fReRPQ2dEBE[2][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
402         fImRPQ2dEBE[2][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
403        }
404       } 
405       // s_{k}(pt,eta) for RP&&POIs // to be improved (clarified)
406       // Remark: final calculation of s_{p,k}(pt,eta) follows after the loop over data bellow:
407       for(Int_t k=0;k<9;k++)
408       {
409        fs2dEBE[2][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k),1.);
410       }
411      } // end of if(fCalculate2DFlow) 
412      */
413       
414     } // end of if(aftsTrack->InPOISelection())
415     
416
417      
418    } // end of if(pTrack->InRPSelection())
419
420   
421   
422    if(aftsTrack->InPOISelection())
423    {
424     dPhi = aftsTrack->Phi();
425     dPt  = aftsTrack->Pt();
426     dEta = aftsTrack->Eta();
427     
428     // 1D (pt)
429     // p_n(m*n,0):   
430     for(Int_t m=0;m<4;m++)
431     {
432      fReRPQ1dEBE[1][0][m][0]->Fill(dPt,TMath::Cos((m+1.)*n*dPhi),1.);
433      fImRPQ1dEBE[1][0][m][0]->Fill(dPt,TMath::Sin((m+1.)*n*dPhi),1.);
434     }
435     // 1D (eta)
436     // p_n(m*n,0):   
437     for(Int_t m=0;m<4;m++)
438     {
439      fReRPQ1dEBE[1][1][m][0]->Fill(dEta,TMath::Cos((m+1.)*n*dPhi),1.);
440      fImRPQ1dEBE[1][1][m][0]->Fill(dEta,TMath::Sin((m+1.)*n*dPhi),1.);
441     }
442     
443     
444     /*
445     // 2D (pt,eta):
446     if(fCalculate2DFlow)
447     {      
448      // p_n(m*n,0):   
449      for(Int_t m=0;m<4;m++)
450      {
451       fReRPQ2dEBE[1][m][0]->Fill(dPt,dEta,TMath::Cos((m+1.)*n*dPhi),1.);
452       fImRPQ2dEBE[1][m][0]->Fill(dPt,dEta,TMath::Sin((m+1.)*n*dPhi),1.);
453      }
454     } // end of if(fCalculate2DFlow)  
455     */
456     
457     
458    } // end of if(pTrack->InPOISelection() )   
459  
460   
461   } else // to if(aftsTrack)
462     {
463      cout<<endl;
464      cout<<" WARNING: no particle! (i.e. aftsTrack is a NULL pointer in AFAWQC::Make().)"<<endl;
465      cout<<endl;       
466     }
467  } // end of for(Int_t i=0;i<nPrim;i++) 
468
469  // calculate the final expressions for S^{M}_{p,k}:
470  for(Int_t p=0;p<8;p++)
471  {
472   for(Int_t k=0;k<9;k++)
473   {
474    (*fSMpk)(p,k)=pow((*fSMpk)(p,k),p+1);
475   }  
476  } 
477  
478  // *****************************
479  // **** CALL THE METHODS *******
480  // *****************************
481  // integrated flow:
482  if(!fEvaluateIntFlowNestedLoops)
483  {
484   if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
485   {
486    if(nRP>1) this->CalculateIntFlowCorrelations(); // without using particle weights
487   } else 
488     {
489      if(nRP>1) this->CalculateIntFlowCorrelationsUsingParticleWeights(); // with using particle weights   
490     } 
491        
492   if(nRP>3) this->CalculateIntFlowProductOfCorrelations();
493   if(nRP>1) this->CalculateIntFlowSumOfEventWeights();
494   if(nRP>1) this->CalculateIntFlowSumOfProductOfEventWeights();
495   if(fApplyCorrectionForNUA && !(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)) // to be improved (enable correction for NUA also when particle weights are used?)
496   {
497    if(nRP>0) this->CalculateIntFlowCorrectionsForNUASinTerms();
498    if(nRP>0) this->CalculateIntFlowCorrectionsForNUACosTerms();
499   }
500  } // end of if(!fEvaluateIntFlowNestedLoops)
501
502  // differential flow:
503  if(!fEvaluateDiffFlowNestedLoops)
504  {
505   if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
506   {
507    if(nRP>1) // to be improved (move this if somewhere else)
508    {
509     // without using particle weights:
510     this->CalculateDiffFlowCorrelations("RP","Pt"); 
511     this->CalculateDiffFlowCorrelations("RP","Eta");
512     this->CalculateDiffFlowCorrelations("POI","Pt");
513     this->CalculateDiffFlowCorrelations("POI","Eta");
514    }  
515   } else
516     {
517      // with using particle weights:   
518      this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Pt"); 
519      this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Eta"); 
520      this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Pt"); 
521      this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Eta"); 
522     } 
523     
524   this->CalculateDiffFlowProductOfCorrelations("RP","Pt");
525   this->CalculateDiffFlowProductOfCorrelations("RP","Eta");
526   this->CalculateDiffFlowProductOfCorrelations("POI","Pt");
527   this->CalculateDiffFlowProductOfCorrelations("POI","Eta");
528   this->CalculateDiffFlowSumOfEventWeights("RP","Pt");
529   this->CalculateDiffFlowSumOfEventWeights("RP","Eta");
530   this->CalculateDiffFlowSumOfEventWeights("POI","Pt");
531   this->CalculateDiffFlowSumOfEventWeights("POI","Eta");
532   this->CalculateDiffFlowSumOfProductOfEventWeights("RP","Pt");
533   this->CalculateDiffFlowSumOfProductOfEventWeights("RP","Eta");
534   this->CalculateDiffFlowSumOfProductOfEventWeights("POI","Pt");
535   this->CalculateDiffFlowSumOfProductOfEventWeights("POI","Eta");
536   if(fApplyCorrectionForNUA && !(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)) // to be improved (enable correction for NUA also when particle weights are used?)
537   {
538    this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Pt");
539    this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Eta");
540    this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Pt");
541    this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Eta");
542    this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Pt");
543    this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Eta");
544    this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Pt");
545    this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Eta");
546   }
547   
548  } // end of if(!fEvaluateDiffFlowNestedLoops)
549
550
551    
552   // with weights:
553   // ... 
554   
555   /*
556   // 2D differential flow
557   if(fCalculate2DFlow)
558   {
559    // without weights:
560    if(nRP>1) this->CalculateCorrelationsForDifferentialFlow2D("RP");
561    if(nRP>1) this->CalculateCorrelationsForDifferentialFlow2D("POI");
562   
563    // with weights:
564    if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
565    {
566     if(nRP>1) this->CalculateWeightedCorrelationsForDifferentialFlow2D("RP");
567     if(nRP>1) this->CalculateWeightedCorrelationsForDifferentialFlow2D("POI");
568    } 
569   } // end of if(fCalculate2DFlow)
570   */
571   
572  
573  // d) Debugging and cross-checking (evaluate nested loops):
574  //  d1) cross-checking results for integrated flow:
575  if(fEvaluateIntFlowNestedLoops)
576  {
577   if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10 
578   {
579    // without using particle weights:
580    if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
581    {
582     // correlations:
583     this->CalculateIntFlowCorrelations(); // from Q-vectors
584     this->EvaluateIntFlowCorrelationsWithNestedLoops(anEvent); // from nested loops (to be improved: do I have to pass here anEvent or not?)
585     // correction for non-uniform acceptance:
586     this->CalculateIntFlowCorrectionsForNUASinTerms(); // from Q-vectors (sin terms)
587     this->CalculateIntFlowCorrectionsForNUACosTerms(); // from Q-vectors (cos terms)
588     this->EvaluateIntFlowCorrectionsForNUAWithNestedLoops(anEvent); // from nested loops (both sin and cos terms)
589    }
590    // using particle weights:
591    if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
592    {
593     // correlations:
594     this->CalculateIntFlowCorrelationsUsingParticleWeights(); // from Q-vectors
595     this->EvaluateIntFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent); // from nested loops (to be improved: do I have to pass here anEvent or not?)
596    }
597   } else if (nPrim>fMaxAllowedMultiplicity) // to if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity)
598     {
599      cout<<endl;
600      cout<<"Skipping the event because multiplicity is "<<nPrim<<". Too high to evaluate nested loops!"<<endl;
601     } else
602       {
603        cout<<endl;
604        cout<<"Skipping the event because multiplicity is "<<nPrim<<"."<<endl;      
605       } 
606  } // end of if(fEvaluateIntFlowNestedLoops) 
607  
608  //  d2) cross-checking results for differential flow:
609  if(fEvaluateDiffFlowNestedLoops)
610  {
611   if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10
612   {
613    // without using particle weights:
614    if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
615    {
616     // reduced correlations:
617     // Q-vectors:
618     this->CalculateDiffFlowCorrelations("RP","Pt");
619     this->CalculateDiffFlowCorrelations("RP","Eta");
620     this->CalculateDiffFlowCorrelations("POI","Pt");
621     this->CalculateDiffFlowCorrelations("POI","Eta");
622     // nested loops:
623     //this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"RP","Pt"); // to be improved (enabled eventually)
624     //this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"RP","Eta"); // to be improved (enabled eventually)
625     this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"POI","Pt"); // to be improved (do I need to pass here anEvent?)
626     this->EvaluateDiffFlowCorrelationsWithNestedLoops(anEvent,"POI","Eta"); // to be improved (do I need to pass here anEvent?)
627     // reduced corrections for non-uniform acceptance:
628     // Q-vectors:
629     this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Pt");
630     this->CalculateDiffFlowCorrectionsForNUASinTerms("RP","Eta");
631     this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Pt");
632     this->CalculateDiffFlowCorrectionsForNUASinTerms("POI","Eta");
633     this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Pt");
634     this->CalculateDiffFlowCorrectionsForNUACosTerms("RP","Eta");
635     this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Pt");
636     this->CalculateDiffFlowCorrectionsForNUACosTerms("POI","Eta");
637     // nested loops:
638     //this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"RP","Pt"); // to be improved (enabled eventually)
639     //this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"RP","Eta"); // to be improved (enabled eventually)
640     this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"POI","Pt"); // to be improved (do I need to pass here anEvent?)
641     this->EvaluateDiffFlowCorrectionTermsForNUAWithNestedLoops(anEvent,"POI","Eta"); // to be improved (do I need to pass here anEvent?)
642    } // end of if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
643    // using particle weights:
644    if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
645    {
646     this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Pt"); 
647     this->CalculateDiffFlowCorrelationsUsingParticleWeights("RP","Eta"); 
648     this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Pt"); 
649     this->CalculateDiffFlowCorrelationsUsingParticleWeights("POI","Eta"); 
650     this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"RP","Pt"); // to be improved (enabled eventually)
651     this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"RP","Eta"); // to be improved (enabled eventually)
652     this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"POI","Pt"); // to be improved (do I need to pass here anEvent?)
653     this->EvaluateDiffFlowCorrelationsWithNestedLoopsUsingParticleWeights(anEvent,"POI","Eta"); // to be improved (do I need to pass here anEvent?)
654    } // end of if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
655   } // end of if(nPrim>0 && nPrim<=fMaxAllowedMultiplicity) // by default fMaxAllowedMultiplicity = 10
656  } // end of if(fEvaluateDiffFlowNestedLoops) 
657  
658  // e) Reset all event by event quantities: 
659  this->ResetEventByEventQuantities();
660  
661 } // end of AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
662
663
664 //================================================================================================================================
665
666
667 void AliFlowAnalysisWithQCumulants::Finish()
668 {
669  // Calculate the final results.
670  //  a) acces the constants;
671  //  b) access the flags;
672  //  c) calculate the final results for integrated flow (without and with weights);
673  //  d) store in AliFlowCommonHistResults and print the final results for integrated flow;
674  //  e) calculate the final results for differential flow (without and with weights);
675  //  f) print the final results for integrated flow obtained from differential flow (to be improved (terminology));
676  //  g) cross-check the results: results from Q-vectors vs results from nested loops
677  
678  // ******************************
679  // **** ACCESS THE CONSTANTS ****
680  // ******************************
681  
682  this->AccessConstants();          
683  
684  if(fCommonHists && fCommonHists->GetHarmonic())
685  {
686   fHarmonic = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1); // to be improved (moved somewhere else)
687  } 
688
689  // **************************
690  // **** ACCESS THE FLAGS ****
691  // **************************    
692  fUsePhiWeights = (Int_t)fUseParticleWeights->GetBinContent(1); 
693  fUsePtWeights = (Int_t)fUseParticleWeights->GetBinContent(2); 
694  fUseEtaWeights = (Int_t)fUseParticleWeights->GetBinContent(3);  
695  fApplyCorrectionForNUA = (Int_t)fIntFlowFlags->GetBinContent(3); 
696  fEvaluateIntFlowNestedLoops = (Int_t)fEvaluateNestedLoops->GetBinContent(1);
697  fEvaluateDiffFlowNestedLoops = (Int_t)fEvaluateNestedLoops->GetBinContent(2); 
698  fCrossCheckInPtBinNo = (Int_t)fEvaluateNestedLoops->GetBinContent(3);
699  fCrossCheckInEtaBinNo = (Int_t)fEvaluateNestedLoops->GetBinContent(4); 
700     
701  // *********************************************************
702  // **** CALCULATE THE FINAL RESULTS FOR INTEGRATED FLOW ****
703  // *********************************************************     
704  
705  this->FinalizeCorrelationsIntFlow();
706  this->CalculateCovariancesIntFlow();
707  this->CalculateCumulantsIntFlow();
708  this->CalculateIntFlow(); 
709
710  if(fApplyCorrectionForNUA && !(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)) // to be improved (reorganized, etc)
711  {
712   this->FinalizeCorrectionTermsForNUAIntFlow();
713   this->CalculateQcumulantsCorrectedForNUAIntFlow();   
714   this->CalculateIntFlowCorrectedForNUA(); 
715  }
716   
717  // ***************************************************************
718  // **** STORE AND PRINT THE FINAL RESULTS FOR INTEGRATED FLOW ****
719  // ***************************************************************
720  
721  this->FillCommonHistResultsIntFlow();  
722   
723  this->PrintFinalResultsForIntegratedFlow("NONAME"); // to be improved (name)
724  
725  // ***********************************************************
726  // **** CALCULATE THE FINAL RESULTS FOR DIFFERENTIAL FLOW ****
727  // ***********************************************************    
728  
729  this->FinalizeReducedCorrelations("RP","Pt"); 
730  this->FinalizeReducedCorrelations("RP","Eta"); 
731  this->FinalizeReducedCorrelations("POI","Pt"); 
732  this->FinalizeReducedCorrelations("POI","Eta");
733  this->CalculateDiffFlowCovariances("RP","Pt");
734  this->CalculateDiffFlowCovariances("RP","Eta");
735  this->CalculateDiffFlowCovariances("POI","Pt");
736  this->CalculateDiffFlowCovariances("POI","Eta");
737  this->CalculateDiffFlowCumulants("RP","Pt");
738  this->CalculateDiffFlowCumulants("RP","Eta");
739  this->CalculateDiffFlowCumulants("POI","Pt");
740  this->CalculateDiffFlowCumulants("POI","Eta");
741  this->CalculateDiffFlow("RP","Pt");
742  this->CalculateDiffFlow("RP","Eta");
743  this->CalculateDiffFlow("POI","Pt");
744  this->CalculateDiffFlow("POI","Eta");
745  
746  if(fApplyCorrectionForNUA && !(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)) // to be improved (reorganized, etc)
747  {
748   this->FinalizeCorrectionTermsForNUADiffFlow("RP","Pt");
749   this->FinalizeCorrectionTermsForNUADiffFlow("RP","Eta");
750   this->FinalizeCorrectionTermsForNUADiffFlow("POI","Pt");
751   this->FinalizeCorrectionTermsForNUADiffFlow("POI","Eta");      
752   this->CalculateDiffFlowCumulantsCorrectedForNUA("RP","Pt");   
753   this->CalculateDiffFlowCumulantsCorrectedForNUA("RP","Eta");   
754   this->CalculateDiffFlowCumulantsCorrectedForNUA("POI","Pt");   
755   this->CalculateDiffFlowCumulantsCorrectedForNUA("POI","Eta");  
756   this->CalculateDiffFlowCorrectedForNUA("RP","Pt"); 
757   this->CalculateDiffFlowCorrectedForNUA("RP","Eta"); 
758   this->CalculateDiffFlowCorrectedForNUA("POI","Pt"); 
759   this->CalculateDiffFlowCorrectedForNUA("POI","Eta"); 
760  }
761  
762  this->CalculateFinalResultsForRPandPOIIntegratedFlow("RP");
763  this->CalculateFinalResultsForRPandPOIIntegratedFlow("POI");
764
765  // *****************************************************************
766  // **** STORE AND PRINT THE FINAL RESULTS FOR DIFFERENTIAL FLOW ****
767  // *****************************************************************
768  this->FillCommonHistResultsDiffFlow("RP");
769  this->FillCommonHistResultsDiffFlow("POI");
770
771  this->PrintFinalResultsForIntegratedFlow("RP"); 
772  this->PrintFinalResultsForIntegratedFlow("POI"); 
773   
774  // g) cross-check the results: results from Q-vectors vs results from nested loops
775  //  g1) integrated flow:
776  if(fEvaluateIntFlowNestedLoops)
777  {
778   if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)) 
779   {
780    this->CrossCheckIntFlowCorrelations();
781    this->CrossCheckIntFlowCorrectionTermsForNUA(); 
782   } else
783     {
784      this->CrossCheckIntFlowCorrelations();     
785      this->CrossCheckIntFlowExtraCorrelations();     
786     }
787  } // end of if(fEvaluateIntFlowNestedLoops)  
788  //  g2) differential flow: 
789  if(fEvaluateDiffFlowNestedLoops) 
790  {
791   // correlations:
792   //this->CrossCheckDiffFlowCorrelations("RP","Pt"); // to be improved (enabled eventually)  
793   //this->CrossCheckDiffFlowCorrelations("RP","Eta"); // to be improved (enabled eventually)  
794   this->CrossCheckDiffFlowCorrelations("POI","Pt");  
795   this->CrossCheckDiffFlowCorrelations("POI","Eta");
796   // correction terms for non-uniform acceptance:
797   //this->CrossCheckDiffFlowCorrectionTermsForNUA("RP","Pt"); // to be improved (enabled eventually)      
798   //this->CrossCheckDiffFlowCorrectionTermsForNUA("RP","Eta"); // to be improved (enabled eventually)      
799   if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)) 
800   {
801    this->CrossCheckDiffFlowCorrectionTermsForNUA("POI","Pt");      
802    this->CrossCheckDiffFlowCorrectionTermsForNUA("POI","Eta");      
803   } 
804  } // end of if(fEvaluateDiffFlowNestedLoops)
805                                                                                                                                                                                                                                                                                                                                    
806 } // end of AliFlowAnalysisWithQCumulants::Finish()
807
808
809 //================================================================================================================================
810
811
812 void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUACosTerms()
813 {
814  // calculate corrections for non-uniform acceptance of the detector for no-name integrated flow (cos terms)
815  
816  // multiplicity:
817  Double_t dMult = (*fSMpk)(0,0);
818  
819  // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
820  Double_t dReQ1n = (*fReQ)(0,0);
821  Double_t dReQ2n = (*fReQ)(1,0);
822  //Double_t dReQ3n = (*fReQ)(2,0);
823  //Double_t dReQ4n = (*fReQ)(3,0);
824  Double_t dImQ1n = (*fImQ)(0,0);
825  Double_t dImQ2n = (*fImQ)(1,0);
826  //Double_t dImQ3n = (*fImQ)(2,0);
827  //Double_t dImQ4n = (*fImQ)(3,0);
828         
829  //                                  *************************************************************
830  //                                  **** corrections for non-uniform acceptance (cos terms): ****
831  //                                  *************************************************************
832  //
833  // Remark 1: corrections for non-uniform acceptance (cos terms) calculated with non-weighted Q-vectors 
834  //           are stored in 1D profile fQCorrectionsCos.
835  // Remark 2: binning of fIntFlowCorrectionTermsForNUAPro[1] is organized as follows:
836  // --------------------------------------------------------------------------------------------------------------------
837  // 1st bin: <<cos(n*(phi1))>> = cosP1n
838  // 2nd bin: <<cos(n*(phi1+phi2))>> = cosP1nP1n
839  // 3rd bin: <<cos(n*(phi1-phi2-phi3))>> = cosP1nM1nM1n
840  // ...
841  // --------------------------------------------------------------------------------------------------------------------
842   
843  // 1-particle:
844  Double_t cosP1n = 0.; // <<cos(n*(phi1))>>
845    
846  if(dMult>0)
847  {
848   cosP1n = dReQ1n/dMult; 
849   
850   // average non-weighted 1-particle correction (cos terms) for non-uniform acceptance for single event:
851   fIntFlowCorrectionTermsForNUAEBE[1]->SetBinContent(1,cosP1n);
852   
853   // final average non-weighted 1-particle correction (cos terms) for non-uniform acceptance for all events:
854   fIntFlowCorrectionTermsForNUAPro[1]->Fill(0.5,cosP1n,dMult);  
855  } 
856  
857  // 2-particle:
858  Double_t cosP1nP1n = 0.; // <<cos(n*(phi1+phi2))>>
859  
860  if(dMult>1)
861  {
862   cosP1nP1n = (pow(dReQ1n,2)-pow(dImQ1n,2)-dReQ2n)/(dMult*(dMult-1)); 
863   
864   // average non-weighted 2-particle correction (cos terms) for non-uniform acceptance for single event:
865   fIntFlowCorrectionTermsForNUAEBE[1]->SetBinContent(2,cosP1nP1n);
866   
867   // final average non-weighted 2-particle correction (cos terms) for non-uniform acceptance for all events:
868   fIntFlowCorrectionTermsForNUAPro[1]->Fill(1.5,cosP1nP1n,dMult*(dMult-1));  
869  } 
870  
871  // 3-particle:
872  Double_t cosP1nM1nM1n = 0.; // <<cos(n*(phi1-phi2-phi3))>>
873  
874  if(dMult>2)
875  {
876   cosP1nM1nM1n = (dReQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))-dReQ1n*dReQ2n-dImQ1n*dImQ2n-2.*(dMult-1)*dReQ1n)
877                / (dMult*(dMult-1)*(dMult-2)); 
878   
879   // average non-weighted 3-particle correction (cos terms) for non-uniform acceptance for single event:
880   fIntFlowCorrectionTermsForNUAEBE[1]->SetBinContent(3,cosP1nM1nM1n);
881   
882   // final average non-weighted 3-particle correction (cos terms) for non-uniform acceptance for all events:
883   fIntFlowCorrectionTermsForNUAPro[1]->Fill(2.5,cosP1nM1nM1n,dMult*(dMult-1)*(dMult-2));  
884  } 
885  
886 } // end of AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUACosTerms()
887
888
889 //================================================================================================================================
890
891
892 void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUASinTerms()
893 {
894  // calculate corrections for non-uniform acceptance of the detector for no-name integrated flow (sin terms)
895  
896  // multiplicity:
897  Double_t dMult = (*fSMpk)(0,0);
898  
899  // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
900  Double_t dReQ1n = (*fReQ)(0,0);
901  Double_t dReQ2n = (*fReQ)(1,0);
902  //Double_t dReQ3n = (*fReQ)(2,0);
903  //Double_t dReQ4n = (*fReQ)(3,0);
904  Double_t dImQ1n = (*fImQ)(0,0);
905  Double_t dImQ2n = (*fImQ)(1,0);
906  //Double_t dImQ3n = (*fImQ)(2,0);
907  //Double_t dImQ4n = (*fImQ)(3,0);
908         
909  //                                  *************************************************************
910  //                                  **** corrections for non-uniform acceptance (sin terms): ****
911  //                                  *************************************************************
912  //
913  // Remark 1: corrections for non-uniform acceptance (sin terms) calculated with non-weighted Q-vectors 
914  //           are stored in 1D profile fQCorrectionsSin.
915  // Remark 2: binning of fIntFlowCorrectionTermsForNUAPro[0] is organized as follows:
916  // --------------------------------------------------------------------------------------------------------------------
917  // 1st bin: <<sin(n*(phi1))>> = sinP1n
918  // 2nd bin: <<sin(n*(phi1+phi2))>> = sinP1nP1n
919  // 3rd bin: <<sin(n*(phi1-phi2-phi3))>> = sinP1nM1nM1n
920  // ...
921  // --------------------------------------------------------------------------------------------------------------------
922  
923  // 1-particle:
924  Double_t sinP1n = 0.; // <sin(n*(phi1))>
925  
926  if(dMult>0)
927  {
928   sinP1n = dImQ1n/dMult; 
929      
930   // average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for single event:
931   fIntFlowCorrectionTermsForNUAEBE[0]->SetBinContent(1,sinP1n);
932   
933   // final average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for all events:   
934   fIntFlowCorrectionTermsForNUAPro[0]->Fill(0.5,sinP1n,dMult);  
935  } 
936  
937  // 2-particle:
938  Double_t sinP1nP1n = 0.; // <<sin(n*(phi1+phi2))>>
939  
940  if(dMult>1)
941  {
942   sinP1nP1n = (2.*dReQ1n*dImQ1n-dImQ2n)/(dMult*(dMult-1)); 
943      
944   // average non-weighted 2-particle correction (sin terms) for non-uniform acceptance for single event:
945   fIntFlowCorrectionTermsForNUAEBE[0]->SetBinContent(2,sinP1nP1n);
946   
947   // final average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for all events:      
948   fIntFlowCorrectionTermsForNUAPro[0]->Fill(1.5,sinP1nP1n,dMult*(dMult-1));  
949  } 
950  
951  // 3-particle:
952  Double_t sinP1nM1nM1n = 0.; // <<sin(n*(phi1-phi2-phi3))>>
953  
954  if(dMult>2)
955  {
956   sinP1nM1nM1n = (-dImQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))+dReQ1n*dImQ2n-dImQ1n*dReQ2n+2.*(dMult-1)*dImQ1n)
957                / (dMult*(dMult-1)*(dMult-2)); 
958   
959   // average non-weighted 3-particle correction (sin terms) for non-uniform acceptance for single event:
960   fIntFlowCorrectionTermsForNUAEBE[0]->SetBinContent(3,sinP1nM1nM1n);
961   
962   // final average non-weighted 3-particle correction (sin terms) for non-uniform acceptance for all events:  
963   fIntFlowCorrectionTermsForNUAPro[0]->Fill(2.5,sinP1nM1nM1n,dMult*(dMult-1)*(dMult-2));  
964  } 
965  
966 } // end of AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrectionsForNUASinTerms()
967
968
969 //================================================================================================================================
970
971
972 void AliFlowAnalysisWithQCumulants::GetOutputHistograms(TList *outputListHistos)
973 {
974  // a) Get pointers for common control and common result histograms and profiles.
975  // b) Get pointers for histograms with particle weights.
976  // c) Get pointers for histograms and profiles relevant for integrated flow.
977  // d) Get pointers for histograms and profiles relevant for differental flow.
978  // e) Get pointers for histograms and profiles holding results obtained with nested loops.
979  
980  if(outputListHistos)
981  {      
982   this->GetPointersForCommonHistograms(outputListHistos); // to be improved (no need to pass here argument, use setter for base list instead)
983   this->GetPointersForParticleWeightsHistograms(outputListHistos); // to be improved (no need to pass here argument, use setter for base list instead)
984   this->GetPointersForIntFlowHistograms(outputListHistos); // to be improved (no need to pass here argument, use setter for base list instead)
985   this->GetPointersForDiffFlowHistograms(outputListHistos); // to be improved (no need to pass here argument, use setter for base list instead)
986   this->GetPointersForNestedLoopsHistograms(outputListHistos); // to be improved (no need to pass here argument, use setter for base list instead)
987  }
988    
989 } // end of void AliFlowAnalysisWithQCumulants::GetOutputHistograms(TList *outputListHistos)
990
991
992 //================================================================================================================================
993
994
995 TProfile* AliFlowAnalysisWithQCumulants::MakePtProjection(TProfile2D *profilePtEta) const
996 {
997  // project 2D profile onto pt axis to get 1D profile
998  
999  Int_t nBinsPt   = profilePtEta->GetNbinsX();
1000  Double_t dPtMin = (profilePtEta->GetXaxis())->GetXmin();
1001  Double_t dPtMax = (profilePtEta->GetXaxis())->GetXmax();
1002  
1003  Int_t nBinsEta   = profilePtEta->GetNbinsY();
1004  
1005  TProfile *profilePt = new TProfile("","",nBinsPt,dPtMin,dPtMax); 
1006  
1007  for(Int_t p=1;p<=nBinsPt;p++)
1008  {
1009   Double_t contentPt = 0.;
1010   Double_t entryPt = 0.;
1011   Double_t spreadPt = 0.;
1012   Double_t sum1 = 0.;
1013   Double_t sum2 = 0.;
1014   Double_t sum3 = 0.;
1015   for(Int_t e=1;e<=nBinsEta;e++)
1016   {
1017    contentPt += (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)))
1018               * (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
1019    entryPt   += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
1020    
1021    sum1 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)))
1022          * (pow(profilePtEta->GetBinError(profilePtEta->GetBin(p,e)),2.)
1023             + pow(profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)),2.)); 
1024    sum2 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
1025    sum3 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)))
1026          * (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)));            
1027   }
1028   if(sum2>0. && sum1/sum2-pow(sum3/sum2,2.) > 0.)
1029   {
1030    spreadPt = pow(sum1/sum2-pow(sum3/sum2,2.),0.5);
1031   }
1032   profilePt->SetBinContent(p,contentPt);
1033   profilePt->SetBinEntries(p,entryPt);
1034   {
1035    profilePt->SetBinError(p,spreadPt);
1036   }
1037   
1038  }
1039  
1040  return profilePt;
1041  
1042 } // end of TProfile* AliFlowAnalysisWithQCumulants::MakePtProjection(TProfile2D *profilePtEta)
1043
1044
1045 //================================================================================================================================
1046
1047
1048 TProfile* AliFlowAnalysisWithQCumulants::MakeEtaProjection(TProfile2D *profilePtEta) const
1049 {
1050  // project 2D profile onto eta axis to get 1D profile
1051  
1052  Int_t nBinsEta   = profilePtEta->GetNbinsY();
1053  Double_t dEtaMin = (profilePtEta->GetYaxis())->GetXmin();
1054  Double_t dEtaMax = (profilePtEta->GetYaxis())->GetXmax();
1055  
1056  Int_t nBinsPt = profilePtEta->GetNbinsX();
1057  
1058  TProfile *profileEta = new TProfile("","",nBinsEta,dEtaMin,dEtaMax); 
1059  
1060  for(Int_t e=1;e<=nBinsEta;e++)
1061  {
1062   Double_t contentEta = 0.;
1063   Double_t entryEta = 0.;
1064   for(Int_t p=1;p<=nBinsPt;p++)
1065   {
1066    contentEta += (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)))
1067               * (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
1068    entryEta   += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
1069   }
1070   profileEta->SetBinContent(e,contentEta);
1071   profileEta->SetBinEntries(e,entryEta);
1072  }
1073  
1074  return profileEta;
1075  
1076 } // end of TProfile* AliFlowAnalysisWithQCumulants::MakeEtaProjection(TProfile2D *profilePtEta)
1077
1078
1079 //================================================================================================================================
1080
1081
1082 void AliFlowAnalysisWithQCumulants::PrintFinalResultsForIntegratedFlow(TString type)
1083 {
1084  // printing on the screen the final results for integrated flow (NONAME, POI and RP) // to be improved (NONAME) 
1085  
1086  Int_t n = fHarmonic; 
1087  
1088  if(type == "NONAME" || type == "RP" || type == "POI")
1089  {
1090   if(!(fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && fCommonHistsResults8th))
1091   {
1092    cout<<"WARNING: fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && fCommonHistsResults8th"<<endl;
1093    cout<<"         is NULL in AFAWQC::PFRFIF() !!!!"<<endl;
1094   }
1095  } else
1096    {
1097     cout<<"WARNING: type in not from {NONAME, RP, POI} in AFAWQC::PFRFIF() !!!!"<<endl;
1098     exit(0);
1099    }
1100  
1101  Double_t dVn[4] = {0.}; // array to hold Vn{2}, Vn{4}, Vn{6} and Vn{8}   
1102  Double_t dVnErr[4] = {0.}; // array to hold errors of Vn{2}, Vn{4}, Vn{6} and Vn{8}   
1103  
1104  if(type == "NONAME")
1105  {
1106   dVn[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinContent(1); 
1107   dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinError(1); 
1108   dVn[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinContent(1); 
1109   dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinError(1); 
1110   dVn[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinContent(1); 
1111   dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinError(1); 
1112   dVn[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinContent(1); 
1113   dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinError(1); 
1114  } else if(type == "RP")
1115    {
1116     dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinContent(1); 
1117     dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinError(1); 
1118     dVn[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinContent(1); 
1119     dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinError(1); 
1120     dVn[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinContent(1); 
1121     dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinError(1); 
1122     dVn[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinContent(1); 
1123     dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinError(1); 
1124    } else if(type == "POI")
1125      {
1126       dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinContent(1); 
1127       dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinError(1); 
1128       dVn[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinContent(1); 
1129       dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinError(1); 
1130       dVn[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinContent(1); 
1131       dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinError(1); 
1132       dVn[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinContent(1); 
1133       dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinError(1); 
1134      }
1135  
1136  TString title = " flow estimates from Q-cumulants"; 
1137  TString subtitle = "    ("; 
1138  
1139  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
1140  {
1141   subtitle.Append(type);
1142   subtitle.Append(", without weights)");
1143  } else  
1144    {
1145     subtitle.Append(type);
1146     subtitle.Append(", with weights)");
1147    }
1148   
1149  cout<<endl;
1150  cout<<"*************************************"<<endl;
1151  cout<<"*************************************"<<endl;
1152  cout<<title.Data()<<endl; 
1153  cout<<subtitle.Data()<<endl; 
1154  cout<<endl;
1155   
1156  for(Int_t i=0;i<4;i++)
1157  {
1158   if(dVn[i]>=0.)
1159   {
1160    cout<<"  v_"<<n<<"{"<<2*(i+1)<<"} = "<<dVn[i]<<" +/- "<<dVnErr[i]<<endl;
1161   }
1162   else
1163   {
1164    cout<<"  v_"<<n<<"{"<<2*(i+1)<<"} = Im"<<endl;
1165   }  
1166  }
1167
1168  cout<<endl;
1169  /*
1170  if(type == "NONAME")
1171  {
1172   cout<<"     nEvts = "<<nEvtsNoName<<", AvM = "<<dMultNoName<<endl; // to be improved
1173  }
1174  else if (type == "RP")
1175  {
1176   cout<<"     nEvts = "<<nEvtsRP<<", AvM = "<<dMultRP<<endl; // to be improved  
1177  } 
1178  else if (type == "POI")
1179  {
1180   cout<<"     nEvts = "<<nEvtsPOI<<", AvM = "<<dMultPOI<<endl; // to be improved  
1181  } 
1182  */
1183  cout<<"*************************************"<<endl;
1184  cout<<"*************************************"<<endl;
1185  cout<<endl; 
1186   
1187 }// end of AliFlowAnalysisWithQCumulants::PrintFinalResultsForIntegratedFlow(TString type="NONAME");
1188
1189
1190 //================================================================================================================================
1191
1192
1193 void AliFlowAnalysisWithQCumulants::WriteHistograms(TString outputFileName)
1194 {
1195  //store the final results in output .root file
1196  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
1197  //output->WriteObject(fHistList, "cobjQC","SingleKey");
1198  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1199  delete output;
1200 }
1201
1202
1203 //================================================================================================================================
1204
1205
1206 void AliFlowAnalysisWithQCumulants::BookCommonHistograms()
1207 {
1208  // Book common control histograms and common histograms for final results.
1209  // common control histogram (ALL events)
1210  TString commonHistsName = "AliFlowCommonHistQC";
1211  commonHistsName += fAnalysisLabel->Data();
1212  fCommonHists = new AliFlowCommonHist(commonHistsName.Data());
1213  fHistList->Add(fCommonHists);  
1214  // common control histogram (for events with 2 and more particles)
1215  TString commonHists2ndOrderName = "AliFlowCommonHist2ndOrderQC";
1216  commonHists2ndOrderName += fAnalysisLabel->Data();
1217  fCommonHists2nd = new AliFlowCommonHist(commonHists2ndOrderName.Data());
1218  fHistList->Add(fCommonHists2nd);  
1219  // common control histogram (for events with 4 and more particles)
1220  TString commonHists4thOrderName = "AliFlowCommonHist4thOrderQC";
1221  commonHists4thOrderName += fAnalysisLabel->Data();
1222  fCommonHists4th = new AliFlowCommonHist(commonHists4thOrderName.Data());
1223  fHistList->Add(fCommonHists4th);  
1224  // common control histogram (for events with 6 and more particles)
1225  TString commonHists6thOrderName = "AliFlowCommonHist6thOrderQC";
1226  commonHists6thOrderName += fAnalysisLabel->Data();
1227  fCommonHists6th = new AliFlowCommonHist(commonHists6thOrderName.Data());
1228  fHistList->Add(fCommonHists6th);  
1229  // common control histogram (for events with 8 and more particles)
1230  TString commonHists8thOrderName = "AliFlowCommonHist8thOrderQC";
1231  commonHists8thOrderName += fAnalysisLabel->Data();
1232  fCommonHists8th = new AliFlowCommonHist(commonHists8thOrderName.Data());
1233  fHistList->Add(fCommonHists8th);    
1234  // common histograms for final results (calculated for events with 2 and more particles)
1235  TString commonHistResults2ndOrderName = "AliFlowCommonHistResults2ndOrderQC";
1236  commonHistResults2ndOrderName += fAnalysisLabel->Data();
1237  fCommonHistsResults2nd = new AliFlowCommonHistResults(commonHistResults2ndOrderName.Data());
1238  fHistList->Add(fCommonHistsResults2nd);  
1239  // common histograms for final results (calculated for events with 4 and more particles)
1240  TString commonHistResults4thOrderName = "AliFlowCommonHistResults4thOrderQC";
1241  commonHistResults4thOrderName += fAnalysisLabel->Data();
1242  fCommonHistsResults4th = new AliFlowCommonHistResults(commonHistResults4thOrderName.Data());
1243  fHistList->Add(fCommonHistsResults4th); 
1244  // common histograms for final results (calculated for events with 6 and more particles)
1245  TString commonHistResults6thOrderName = "AliFlowCommonHistResults6thOrderQC";
1246  commonHistResults6thOrderName += fAnalysisLabel->Data();
1247  fCommonHistsResults6th = new AliFlowCommonHistResults(commonHistResults6thOrderName.Data());
1248  fHistList->Add(fCommonHistsResults6th);  
1249  // common histograms for final results (calculated for events with 8 and more particles)
1250  TString commonHistResults8thOrderName = "AliFlowCommonHistResults8thOrderQC";
1251  commonHistResults8thOrderName += fAnalysisLabel->Data();
1252  fCommonHistsResults8th = new AliFlowCommonHistResults(commonHistResults8thOrderName.Data());
1253  fHistList->Add(fCommonHistsResults8th); 
1254  
1255 } // end of void AliFlowAnalysisWithQCumulants::BookCommonHistograms()
1256
1257
1258 //================================================================================================================================
1259
1260
1261 void AliFlowAnalysisWithQCumulants::BookAndFillWeightsHistograms()
1262 {
1263  // book and fill histograms which hold phi, pt and eta weights
1264
1265  if(!fWeightsList)
1266  {
1267   cout<<"WARNING: fWeightsList is NULL in AFAWQC::BAFWH() !!!!"<<endl;
1268   exit(0);  
1269  }
1270     
1271  TString fUseParticleWeightsName = "fUseParticleWeightsQC";
1272  fUseParticleWeightsName += fAnalysisLabel->Data();
1273  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
1274  fUseParticleWeights->SetLabelSize(0.06);
1275  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
1276  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
1277  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
1278  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
1279  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
1280  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
1281  fWeightsList->Add(fUseParticleWeights); 
1282   
1283  if(fUsePhiWeights)
1284  {
1285   if(fWeightsList->FindObject("phi_weights"))
1286   {
1287    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
1288    if(fPhiWeights->GetBinWidth(1) != fPhiBinWidth)
1289    {
1290     cout<<"WARNING: fPhiWeights->GetBinWidth(1) != fPhiBinWidth in AFAWQC::BAFWH() !!!!        "<<endl;
1291     cout<<"         This indicates inconsistent binning in phi histograms throughout the code."<<endl;
1292     exit(0);
1293    }
1294   } else 
1295     {
1296      cout<<"WARNING: fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
1297      exit(0);
1298     }
1299  } // end of if(fUsePhiWeights)
1300  
1301  if(fUsePtWeights) 
1302  {
1303   if(fWeightsList->FindObject("pt_weights"))
1304   {
1305    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
1306    if(fPtWeights->GetBinWidth(1) != fPtBinWidth)
1307    {
1308     cout<<"WARNING: fPtWeights->GetBinWidth(1) != fPtBinWidth in AFAWQC::BAFWH() !!!!         "<<endl;
1309     cout<<"         This indicates insconsistent binning in pt histograms throughout the code."<<endl;
1310     exit(0);
1311    }
1312   } else 
1313     {
1314      cout<<"WARNING: fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
1315      exit(0);
1316     }
1317  } // end of if(fUsePtWeights)    
1318
1319  if(fUseEtaWeights) 
1320  {
1321   if(fWeightsList->FindObject("eta_weights"))
1322   {
1323    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
1324    if(fEtaWeights->GetBinWidth(1) != fEtaBinWidth)
1325    {
1326     cout<<"WARNING: fEtaWeights->GetBinWidth(1) != fEtaBinWidth in AFAWQC::BAFWH() !!!!        "<<endl;
1327     cout<<"         This indicates insconsistent binning in eta histograms throughout the code."<<endl;
1328     exit(0);
1329    }
1330   } else 
1331     {
1332      cout<<"WARNING: fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
1333      exit(0);
1334     }
1335  } // end of if(fUseEtaWeights)
1336  
1337 } // end of AliFlowAnalysisWithQCumulants::BookAndFillWeightsHistograms()
1338
1339
1340 //================================================================================================================================
1341
1342
1343 void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
1344 {
1345  // Book all objects for integrated flow:
1346  //  a) Book profile to hold all flags for integrated flow.
1347  //  b) Book event-by-event quantities.
1348  //  c) Book profiles. // to be improved (comment)
1349  //  d) Book histograms holding the final results.
1350  
1351  TString sinCosFlag[2] = {"sin","cos"}; // to be improved (should I promote this to data members?)
1352  TString powerFlag[2] = {"linear","quadratic"}; // to be improved (should I promote this to data members?)
1353  
1354  // a) Book profile to hold all flags for integrated flow:
1355  TString intFlowFlagsName = "fIntFlowFlags";
1356  intFlowFlagsName += fAnalysisLabel->Data();
1357  fIntFlowFlags = new TProfile(intFlowFlagsName.Data(),"Flags for Integrated Flow",3,0,3);
1358  fIntFlowFlags->SetTickLength(-0.01,"Y");
1359  fIntFlowFlags->SetMarkerStyle(25);
1360  fIntFlowFlags->SetLabelSize(0.05);
1361  fIntFlowFlags->SetLabelOffset(0.02,"Y");
1362  (fIntFlowFlags->GetXaxis())->SetBinLabel(1,"Particle Weights");
1363  (fIntFlowFlags->GetXaxis())->SetBinLabel(2,"Event Weights");
1364  (fIntFlowFlags->GetXaxis())->SetBinLabel(3,"Corrected for NUA?");
1365  fIntFlowList->Add(fIntFlowFlags);
1366
1367  // b) Book event-by-event quantities:
1368  // Re[Q_{m*n,k}], Im[Q_{m*n,k}] and S_{p,k}^M: 
1369  fReQ  = new TMatrixD(4,9);
1370  fImQ  = new TMatrixD(4,9);
1371  fSMpk = new TMatrixD(8,9);
1372  // average correlations <2>, <4>, <6> and <8> for single event (bining is the same as in fIntFlowCorrelationsPro and fIntFlowCorrelationsHist):
1373  TString intFlowCorrelationsEBEName = "fIntFlowCorrelationsEBE";
1374  intFlowCorrelationsEBEName += fAnalysisLabel->Data();
1375  fIntFlowCorrelationsEBE = new TH1D(intFlowCorrelationsEBEName.Data(),intFlowCorrelationsEBEName.Data(),4,0,4);
1376  // weights for average correlations <2>, <4>, <6> and <8> for single event:
1377  TString intFlowEventWeightsForCorrelationsEBEName = "fIntFlowEventWeightsForCorrelationsEBE";
1378  intFlowEventWeightsForCorrelationsEBEName += fAnalysisLabel->Data();
1379  fIntFlowEventWeightsForCorrelationsEBE = new TH1D(intFlowEventWeightsForCorrelationsEBEName.Data(),intFlowEventWeightsForCorrelationsEBEName.Data(),4,0,4);
1380  // average all correlations for single event (bining is the same as in fIntFlowCorrelationsAllPro and fIntFlowCorrelationsAllHist):
1381  TString intFlowCorrelationsAllEBEName = "fIntFlowCorrelationsAllEBE";
1382  intFlowCorrelationsAllEBEName += fAnalysisLabel->Data();
1383  fIntFlowCorrelationsAllEBE = new TH1D(intFlowCorrelationsAllEBEName.Data(),intFlowCorrelationsAllEBEName.Data(),32,0,32);
1384  // average correction terms for non-uniform acceptance for single event 
1385  // (binning is the same as in fIntFlowCorrectionTermsForNUAPro[2] and fIntFlowCorrectionTermsForNUAHist[2]):
1386  TString fIntFlowCorrectionTermsForNUAEBEName = "fIntFlowCorrectionTermsForNUAEBE";
1387  fIntFlowCorrectionTermsForNUAEBEName += fAnalysisLabel->Data();
1388  for(Int_t sc=0;sc<2;sc++) // sin or cos terms
1389  {
1390   fIntFlowCorrectionTermsForNUAEBE[sc] = new TH1D(Form("%s: %s terms",fIntFlowCorrectionTermsForNUAEBEName.Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),10,0,10);  
1391  }
1392  
1393  // c) Book profiles: // to be improved (comment)
1394  // profile to hold average multiplicities and number of events for events with nRP>=0, nRP>=1, ... , and nRP>=8:
1395  TString avMultiplicityName = "fAvMultiplicity";
1396  avMultiplicityName += fAnalysisLabel->Data();
1397  fAvMultiplicity = new TProfile(avMultiplicityName.Data(),"Average Multiplicities of RPs",9,0,9);
1398  fAvMultiplicity->SetTickLength(-0.01,"Y");
1399  fAvMultiplicity->SetMarkerStyle(25);
1400  fAvMultiplicity->SetLabelSize(0.05);
1401  fAvMultiplicity->SetLabelOffset(0.02,"Y");
1402  fAvMultiplicity->SetYTitle("Average Multiplicity");
1403  (fAvMultiplicity->GetXaxis())->SetBinLabel(1,"all evts");
1404  (fAvMultiplicity->GetXaxis())->SetBinLabel(2,"n_{RP} #geq 1");
1405  (fAvMultiplicity->GetXaxis())->SetBinLabel(3,"n_{RP} #geq 2");
1406  (fAvMultiplicity->GetXaxis())->SetBinLabel(4,"n_{RP} #geq 3");
1407  (fAvMultiplicity->GetXaxis())->SetBinLabel(5,"n_{RP} #geq 4");
1408  (fAvMultiplicity->GetXaxis())->SetBinLabel(6,"n_{RP} #geq 5");
1409  (fAvMultiplicity->GetXaxis())->SetBinLabel(7,"n_{RP} #geq 6");
1410  (fAvMultiplicity->GetXaxis())->SetBinLabel(8,"n_{RP} #geq 7");
1411  (fAvMultiplicity->GetXaxis())->SetBinLabel(9,"n_{RP} #geq 8");
1412  fIntFlowProfiles->Add(fAvMultiplicity);
1413  // average correlations <<2>>, <<4>>, <<6>> and <<8>> for all events (with wrong errors!):
1414  TString intFlowCorrelationsProName = "fIntFlowCorrelationsPro";
1415  intFlowCorrelationsProName += fAnalysisLabel->Data();
1416  fIntFlowCorrelationsPro = new TProfile(intFlowCorrelationsProName.Data(),"Average correlations for all events",4,0,4,"s");
1417  fIntFlowCorrelationsPro->SetTickLength(-0.01,"Y");
1418  fIntFlowCorrelationsPro->SetMarkerStyle(25);
1419  fIntFlowCorrelationsPro->SetLabelSize(0.06);
1420  fIntFlowCorrelationsPro->SetLabelOffset(0.01,"Y");
1421  (fIntFlowCorrelationsPro->GetXaxis())->SetBinLabel(1,"<<2>>");
1422  (fIntFlowCorrelationsPro->GetXaxis())->SetBinLabel(2,"<<4>>");
1423  (fIntFlowCorrelationsPro->GetXaxis())->SetBinLabel(3,"<<6>>");
1424  (fIntFlowCorrelationsPro->GetXaxis())->SetBinLabel(4,"<<8>>");
1425  fIntFlowProfiles->Add(fIntFlowCorrelationsPro);
1426  // averaged all correlations for all events (with wrong errors!):
1427  TString intFlowCorrelationsAllProName = "fIntFlowCorrelationsAllPro";
1428  intFlowCorrelationsAllProName += fAnalysisLabel->Data();
1429  fIntFlowCorrelationsAllPro = new TProfile(intFlowCorrelationsAllProName.Data(),"Average correlations for all events",32,0,32,"s");
1430  fIntFlowCorrelationsAllPro->SetTickLength(-0.01,"Y");
1431  fIntFlowCorrelationsAllPro->SetMarkerStyle(25);
1432  fIntFlowCorrelationsAllPro->SetLabelSize(0.03);
1433  fIntFlowCorrelationsAllPro->SetLabelOffset(0.01,"Y");
1434  // 2-p correlations:
1435  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
1436  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
1437  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
1438  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
1439  // 3-p correlations:
1440  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
1441  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
1442  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
1443  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
1444  // 4-p correlations:
1445  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
1446  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
1447  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
1448  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
1449  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
1450  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
1451  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
1452  // 5-p correlations:
1453  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
1454  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
1455  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
1456  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
1457  // 6-p correlations:
1458  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
1459  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
1460  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
1461  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
1462  // 7-p correlations:  
1463  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
1464  // 8-p correlations:
1465  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
1466  fIntFlowProfiles->Add(fIntFlowCorrelationsAllPro);
1467  // when particle weights are used some extra correlations appear:
1468  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights) 
1469  {
1470   TString intFlowExtraCorrelationsProName = "fIntFlowExtraCorrelationsPro";
1471   intFlowExtraCorrelationsProName += fAnalysisLabel->Data();
1472   fIntFlowExtraCorrelationsPro = new TProfile(intFlowExtraCorrelationsProName.Data(),"Average extra correlations for all events",100,0,100,"s");
1473   fIntFlowExtraCorrelationsPro->SetTickLength(-0.01,"Y");
1474   fIntFlowExtraCorrelationsPro->SetMarkerStyle(25);
1475   fIntFlowExtraCorrelationsPro->SetLabelSize(0.03);
1476   fIntFlowExtraCorrelationsPro->SetLabelOffset(0.01,"Y");
1477   // extra 2-p correlations:
1478   (fIntFlowExtraCorrelationsPro->GetXaxis())->SetBinLabel(1,"<<w1^3 w2 cos(n*(phi1-phi2))>>");
1479   (fIntFlowExtraCorrelationsPro->GetXaxis())->SetBinLabel(2,"<<w1 w2 w3^2 cos(n*(phi1-phi2))>>");
1480   // ...
1481   fIntFlowProfiles->Add(fIntFlowExtraCorrelationsPro);
1482  } // end of if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
1483  // average product of correlations <2>, <4>, <6> and <8>:  
1484  TString intFlowProductOfCorrelationsProName = "fIntFlowProductOfCorrelationsPro";
1485  intFlowProductOfCorrelationsProName += fAnalysisLabel->Data();
1486  fIntFlowProductOfCorrelationsPro = new TProfile(intFlowProductOfCorrelationsProName.Data(),"Average products of correlations",6,0,6);
1487  fIntFlowProductOfCorrelationsPro->SetTickLength(-0.01,"Y");
1488  fIntFlowProductOfCorrelationsPro->SetMarkerStyle(25); 
1489  fIntFlowProductOfCorrelationsPro->SetLabelSize(0.05);
1490  fIntFlowProductOfCorrelationsPro->SetLabelOffset(0.01,"Y");
1491  (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(1,"<<2><4>>");
1492  (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(2,"<<2><6>>");
1493  (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(3,"<<2><8>>");
1494  (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(4,"<<4><6>>");
1495  (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(5,"<<4><8>>");
1496  (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(6,"<<6><8>>");
1497  fIntFlowProfiles->Add(fIntFlowProductOfCorrelationsPro);
1498  // average correction terms for non-uniform acceptance (with wrong errors!):
1499  for(Int_t sc=0;sc<2;sc++) // sin or cos terms
1500  {
1501   TString intFlowCorrectionTermsForNUAProName = "fIntFlowCorrectionTermsForNUAPro";
1502   intFlowCorrectionTermsForNUAProName += fAnalysisLabel->Data();
1503   fIntFlowCorrectionTermsForNUAPro[sc] = new TProfile(Form("%s: %s terms",intFlowCorrectionTermsForNUAProName.Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),10,0,10,"s");
1504   fIntFlowCorrectionTermsForNUAPro[sc]->SetTickLength(-0.01,"Y");
1505   fIntFlowCorrectionTermsForNUAPro[sc]->SetMarkerStyle(25);
1506   fIntFlowCorrectionTermsForNUAPro[sc]->SetLabelSize(0.03);
1507   fIntFlowCorrectionTermsForNUAPro[sc]->SetLabelOffset(0.01,"Y");
1508   // 1-particle terms:
1509   (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(1,Form("<<%s(n(phi1))>>",sinCosFlag[sc].Data()));
1510   // 2-particle terms:
1511   (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(2,Form("<<%s(n(phi1+phi2))>>",sinCosFlag[sc].Data()));  
1512   // 3-particle terms:
1513   (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(3,Form("<<%s(n(phi1-phi2-phi3))>>",sinCosFlag[sc].Data()));  
1514   // ... 
1515   fIntFlowProfiles->Add(fIntFlowCorrectionTermsForNUAPro[sc]);
1516  } // end of for(Int_t sc=0;sc<2;sc++) 
1517  
1518  // d) Book histograms holding the final results:
1519  // average correlations <<2>>, <<4>>, <<6>> and <<8>> for all events (with correct errors!):
1520  TString intFlowCorrelationsHistName = "fIntFlowCorrelationsHist";
1521  intFlowCorrelationsHistName += fAnalysisLabel->Data();
1522  fIntFlowCorrelationsHist = new TH1D(intFlowCorrelationsHistName.Data(),"Average correlations for all events",4,0,4);
1523  fIntFlowCorrelationsHist->SetTickLength(-0.01,"Y");
1524  fIntFlowCorrelationsHist->SetMarkerStyle(25);
1525  fIntFlowCorrelationsHist->SetLabelSize(0.06);
1526  fIntFlowCorrelationsHist->SetLabelOffset(0.01,"Y");
1527  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(1,"<<2>>");
1528  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(2,"<<4>>");
1529  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(3,"<<6>>");
1530  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(4,"<<8>>");
1531  fIntFlowResults->Add(fIntFlowCorrelationsHist);
1532  // average all correlations for all events (with correct errors!):
1533  TString intFlowCorrelationsAllHistName = "fIntFlowCorrelationsAllHist";
1534  intFlowCorrelationsAllHistName += fAnalysisLabel->Data();
1535  fIntFlowCorrelationsAllHist = new TH1D(intFlowCorrelationsAllHistName.Data(),"Average correlations for all events",32,0,32);
1536  fIntFlowCorrelationsAllHist->SetTickLength(-0.01,"Y");
1537  fIntFlowCorrelationsAllHist->SetMarkerStyle(25);
1538  fIntFlowCorrelationsAllHist->SetLabelSize(0.03);
1539  fIntFlowCorrelationsAllHist->SetLabelOffset(0.01,"Y");
1540  // 2-p correlations:
1541  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
1542  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
1543  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
1544  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
1545  // 3-p correlations:
1546  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
1547  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
1548  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
1549  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
1550  // 4-p correlations:
1551  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
1552  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
1553  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
1554  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
1555  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
1556  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
1557  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
1558  // 5-p correlations:
1559  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
1560  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
1561  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
1562  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
1563  // 6-p correlations:
1564  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
1565  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
1566  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
1567  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
1568  // 7-p correlations:  
1569  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
1570  // 8-p correlations:
1571  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
1572  fIntFlowResults->Add(fIntFlowCorrelationsAllHist);
1573  // average correction terms for non-uniform acceptance (with correct errors!):
1574  for(Int_t sc=0;sc<2;sc++) // sin or cos terms
1575  {
1576   TString intFlowCorrectionTermsForNUAHistName = "fIntFlowCorrectionTermsForNUAHist";
1577   intFlowCorrectionTermsForNUAHistName += fAnalysisLabel->Data();
1578   fIntFlowCorrectionTermsForNUAHist[sc] = new TH1D(Form("%s: %s terms",intFlowCorrectionTermsForNUAHistName.Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),10,0,10);
1579   fIntFlowCorrectionTermsForNUAHist[sc]->SetTickLength(-0.01,"Y");
1580   fIntFlowCorrectionTermsForNUAHist[sc]->SetMarkerStyle(25);
1581   fIntFlowCorrectionTermsForNUAHist[sc]->SetLabelSize(0.03);
1582   fIntFlowCorrectionTermsForNUAHist[sc]->SetLabelOffset(0.01,"Y");
1583   // ......................................................................... 
1584   // 1-p terms:
1585   (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(1,Form("%s(n(#phi_{1}))>",sinCosFlag[sc].Data()));
1586   // 2-p terms:
1587   // 3-p terms:
1588   // ...
1589   // ......................................................................... 
1590   fIntFlowResults->Add(fIntFlowCorrectionTermsForNUAHist[sc]);
1591  } // end of for(Int_t sc=0;sc<2;sc++) 
1592  // covariances (multiplied with weight dependent prefactor):
1593  TString intFlowCovariancesName = "fIntFlowCovariances";
1594  intFlowCovariancesName += fAnalysisLabel->Data();
1595  fIntFlowCovariances = new TH1D(intFlowCovariancesName.Data(),"Covariances (multiplied with weight dependent prefactor)",6,0,6);
1596  fIntFlowCovariances->SetLabelSize(0.04);
1597  fIntFlowCovariances->SetMarkerStyle(25);
1598  (fIntFlowCovariances->GetXaxis())->SetBinLabel(1,"Cov(<2>,<4>)");
1599  (fIntFlowCovariances->GetXaxis())->SetBinLabel(2,"Cov(<2>,<6>)");
1600  (fIntFlowCovariances->GetXaxis())->SetBinLabel(3,"Cov(<2>,<8>)");
1601  (fIntFlowCovariances->GetXaxis())->SetBinLabel(4,"Cov(<4>,<6>)");
1602  (fIntFlowCovariances->GetXaxis())->SetBinLabel(5,"Cov(<4>,<8>)");
1603  (fIntFlowCovariances->GetXaxis())->SetBinLabel(6,"Cov(<6>,<8>)");  
1604  fIntFlowResults->Add(fIntFlowCovariances);
1605  // sum of linear and quadratic event weights for <2>, <4>, <6> and <8>:
1606  TString intFlowSumOfEventWeightsName = "fIntFlowSumOfEventWeights";
1607  intFlowSumOfEventWeightsName += fAnalysisLabel->Data();
1608  for(Int_t power=0;power<2;power++)
1609  {
1610   fIntFlowSumOfEventWeights[power] = new TH1D(Form("%s: %s",intFlowSumOfEventWeightsName.Data(),powerFlag[power].Data()),Form("Sum of %s event weights for correlations",powerFlag[power].Data()),4,0,4);
1611   fIntFlowSumOfEventWeights[power]->SetLabelSize(0.05);
1612   fIntFlowSumOfEventWeights[power]->SetMarkerStyle(25);
1613   if(power == 0)
1614   {
1615    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{<2>}");
1616    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{<4>}");
1617    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{<6>}");
1618    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{<8>}");
1619   } else if (power == 1) 
1620     {
1621      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{<2>}^{2}");
1622      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{<4>}^{2}");
1623      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{<6>}^{2}");
1624      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{<8>}^{2}");
1625     }
1626   fIntFlowResults->Add(fIntFlowSumOfEventWeights[power]);
1627  } 
1628  // sum of products of event weights for correlations <2>, <4>, <6> and <8>:  
1629  TString intFlowSumOfProductOfEventWeightsName = "fIntFlowSumOfProductOfEventWeights";
1630  intFlowSumOfProductOfEventWeightsName += fAnalysisLabel->Data();
1631  fIntFlowSumOfProductOfEventWeights = new TH1D(intFlowSumOfProductOfEventWeightsName.Data(),"Sum of product of event weights for correlations",6,0,6);
1632  fIntFlowSumOfProductOfEventWeights->SetLabelSize(0.05);
1633  fIntFlowSumOfProductOfEventWeights->SetMarkerStyle(25);
1634  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{<2>} w_{<4>}");
1635  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{<2>} w_{<6>}");
1636  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{<2>} w_{<8>}");
1637  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{<4>} w_{<6>}");
1638  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(5,"#sum_{i=1}^{N} w_{<4>} w_{<8>}");
1639  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(6,"#sum_{i=1}^{N} w_{<6>} w_{<8>}");
1640  fIntFlowResults->Add(fIntFlowSumOfProductOfEventWeights);
1641  // final results for integrated Q-cumulants:
1642  TString intFlowQcumulantsName = "fIntFlowQcumulants";
1643  intFlowQcumulantsName += fAnalysisLabel->Data();
1644  fIntFlowQcumulants = new TH1D(intFlowQcumulantsName.Data(),"Integrated Q-cumulants",4,0,4);
1645  fIntFlowQcumulants->SetLabelSize(0.05);
1646  fIntFlowQcumulants->SetMarkerStyle(25);
1647  (fIntFlowQcumulants->GetXaxis())->SetBinLabel(1,"QC{2}");
1648  (fIntFlowQcumulants->GetXaxis())->SetBinLabel(2,"QC{4}");
1649  (fIntFlowQcumulants->GetXaxis())->SetBinLabel(3,"QC{6}");
1650  (fIntFlowQcumulants->GetXaxis())->SetBinLabel(4,"QC{8}");
1651  fIntFlowResults->Add(fIntFlowQcumulants);
1652  // final integrated flow estimates from Q-cumulants:
1653  TString intFlowName = "fIntFlow";
1654  intFlowName += fAnalysisLabel->Data();  
1655  // integrated flow from Q-cumulants:
1656  fIntFlow = new TH1D(intFlowName.Data(),"Integrated flow estimates from Q-cumulants",4,0,4);
1657  fIntFlow->SetLabelSize(0.05);
1658  fIntFlow->SetMarkerStyle(25);
1659  (fIntFlow->GetXaxis())->SetBinLabel(1,"v_{2}{2,QC}");
1660  (fIntFlow->GetXaxis())->SetBinLabel(2,"v_{2}{4,QC}");
1661  (fIntFlow->GetXaxis())->SetBinLabel(3,"v_{2}{6,QC}");
1662  (fIntFlow->GetXaxis())->SetBinLabel(4,"v_{2}{8,QC}");
1663  fIntFlowResults->Add(fIntFlow);
1664
1665  /* // to be improved (removed):
1666   // final average weighted multi-particle correlations for all events calculated from Q-vectors
1667   fQCorrelations[1] = new TProfile("Weighted correlations","final average multi-particle correlations from weighted Q-vectors",200,0,200,"s");
1668   fQCorrelations[1]->SetTickLength(-0.01,"Y");
1669   fQCorrelations[1]->SetMarkerStyle(25);
1670   fQCorrelations[1]->SetLabelSize(0.03);
1671   fQCorrelations[1]->SetLabelOffset(0.01,"Y");
1672   // 2-particle correlations:
1673   (fQCorrelations[1]->GetXaxis())->SetBinLabel(1,"<w_{1}w_{2}cos(n(#phi_{1}-#phi_{2}))>");
1674   (fQCorrelations[1]->GetXaxis())->SetBinLabel(2,"<w_{1}^{2}w_{2}^{2}cos(2n(#phi_{1}-#phi_{2}))>");
1675   (fQCorrelations[1]->GetXaxis())->SetBinLabel(3,"<w_{1}^{3}w_{2}^{3}cos(3n(#phi_{1}-#phi_{2}))>");
1676   (fQCorrelations[1]->GetXaxis())->SetBinLabel(4,"<w_{1}^{4}w_{2}^{4}cos(4n(#phi_{1}-#phi_{2}))>");
1677   (fQCorrelations[1]->GetXaxis())->SetBinLabel(5,"<w_{1}^{3}w_{2}cos(n(#phi_{1}-#phi_{2}))>");
1678   (fQCorrelations[1]->GetXaxis())->SetBinLabel(6,"<w_{1}^{2}w_{2}w_{3}cos(n(#phi_{1}-#phi_{2}))>");
1679   // 3-particle correlations:
1680   (fQCorrelations[1]->GetXaxis())->SetBinLabel(21,"<w_{1}w_{2}w_{3}^{2}cos(n(2#phi_{1}-#phi_{2}-#phi_{3}))>");
1681   // 4-particle correlations:
1682   (fQCorrelations[1]->GetXaxis())->SetBinLabel(41,"<w_{1}w_{2}w_{3}w_{4}cos(n(#phi_{1}+#phi_{2}-#phi_{3}-#phi_{4}))>");
1683   // add fQCorrelations[1] to the list fIntFlowList:
1684   fIntFlowList->Add(fQCorrelations[1]); 
1685  */
1686   
1687 } // end of AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
1688
1689
1690 //================================================================================================================================
1691
1692
1693 void AliFlowAnalysisWithQCumulants::InitializeArraysForNestedLoops()
1694 {
1695  // Initialize arrays of all objects relevant for calculations with nested loops.
1696  
1697  // integrated flow:
1698  for(Int_t sc=0;sc<2;sc++) // sin or cos terms
1699  {
1700   fIntFlowDirectCorrectionTermsForNUA[sc] = NULL;
1701  } 
1702
1703  // differential flow:  
1704  // correlations:
1705  for(Int_t t=0;t<2;t++) // type: RP or POI
1706  { 
1707   for(Int_t pe=0;pe<2;pe++) // pt or eta
1708   {
1709    for(Int_t ci=0;ci<4;ci++) // correlation index
1710    {
1711     fDiffFlowDirectCorrelations[t][pe][ci] = NULL;
1712    } // end of for(Int_t ci=0;ci<4;ci++) // correlation index  
1713   } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
1714  } // end of for(Int_t t=0;t<2;t++) // type: RP or POI
1715  // correction terms for non-uniform acceptance:
1716  for(Int_t t=0;t<2;t++) // type: RP or POI
1717  { 
1718   for(Int_t pe=0;pe<2;pe++) // pt or eta
1719   {
1720    for(Int_t sc=0;sc<2;sc++) // sin or cos terms
1721    {
1722     for(Int_t cti=0;cti<9;cti++) // correction term index
1723     {
1724      fDiffFlowDirectCorrectionTermsForNUA[t][pe][sc][cti] = NULL;
1725     }   
1726    }
1727   } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
1728  } // end of for(Int_t t=0;t<2;t++) // type: RP or POI
1729
1730
1731 } // end of void AliFlowAnalysisWithQCumulants::InitializeArraysForNestedLoops()
1732
1733
1734 //================================================================================================================================
1735
1736
1737 void AliFlowAnalysisWithQCumulants::BookEverythingForNestedLoops()
1738 {
1739  // Book all objects relevant for calculations with nested loops.
1740  
1741  TString sinCosFlag[2] = {"sin","cos"}; // to be improved (should I promote this to data members?)
1742  TString typeFlag[2] = {"RP","POI"}; // to be improved (should I promote this to data members?)
1743  TString ptEtaFlag[2] = {"p_{T}","#eta"}; // to be improved (should I promote this to data members?)
1744  TString reducedCorrelationIndex[4] = {"<2'>","<4'>","<6'>","<8'>"}; // to be improved (should I promote this to data members?)
1745  Double_t lowerPtEtaEdge[2] = {fPtMin+(fCrossCheckInPtBinNo-1)*fPtBinWidth,fEtaMin+(fCrossCheckInEtaBinNo-1)*fEtaBinWidth};
1746  Double_t upperPtEtaEdge[2] = {fPtMin+fCrossCheckInPtBinNo*fPtBinWidth,fEtaMin+fCrossCheckInEtaBinNo*fEtaBinWidth};
1747
1748  TString evaluateNestedLoopsName = "fEvaluateNestedLoops";
1749  evaluateNestedLoopsName += fAnalysisLabel->Data();
1750  fEvaluateNestedLoops = new TProfile(evaluateNestedLoopsName.Data(),"Flags for nested loops",4,0,4);
1751  fEvaluateNestedLoops->SetLabelSize(0.03);
1752  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(1,"fEvaluateIntFlowNestedLoops");
1753  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(2,"fEvaluateDiffFlowNestedLoops");
1754  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(3,"fCrossCheckInPtBinNo");
1755  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(4,"fCrossCheckInEtaBinNo");
1756  fEvaluateNestedLoops->Fill(0.5,(Int_t)fEvaluateIntFlowNestedLoops);
1757  fEvaluateNestedLoops->Fill(1.5,(Int_t)fEvaluateDiffFlowNestedLoops);
1758  fEvaluateNestedLoops->Fill(2.5,fCrossCheckInPtBinNo);
1759  fEvaluateNestedLoops->Fill(3.5,fCrossCheckInEtaBinNo);
1760  fNestedLoopsList->Add(fEvaluateNestedLoops);
1761  // nested loops for integrated flow:
1762  if(fEvaluateIntFlowNestedLoops)
1763  {
1764   // correlations:
1765   TString intFlowDirectCorrelationsName = "fIntFlowDirectCorrelations";
1766   intFlowDirectCorrelationsName += fAnalysisLabel->Data();
1767   fIntFlowDirectCorrelations = new TProfile(intFlowDirectCorrelationsName.Data(),"Multiparticle correlations calculated with nested loops (for int. flow)",32,0,32,"s");
1768   fNestedLoopsList->Add(fIntFlowDirectCorrelations);
1769   if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
1770   {
1771    TString intFlowExtraDirectCorrelationsName = "fIntFlowExtraDirectCorrelations";
1772    intFlowExtraDirectCorrelationsName += fAnalysisLabel->Data();
1773    fIntFlowExtraDirectCorrelations = new TProfile(intFlowExtraDirectCorrelationsName.Data(),"Extra multiparticle correlations calculated with nested loops (for int. flow)",100,0,100,"s");
1774    fNestedLoopsList->Add(fIntFlowExtraDirectCorrelations);  
1775   } // end of if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
1776   // correction terms for non-uniform acceptance:
1777   for(Int_t sc=0;sc<2;sc++) // sin or cos terms
1778   {
1779    TString intFlowDirectCorrectionTermsForNUAName = "fIntFlowDirectCorrectionTermsForNUA";
1780    intFlowDirectCorrectionTermsForNUAName += fAnalysisLabel->Data();
1781    fIntFlowDirectCorrectionTermsForNUA[sc] = new TProfile(Form("%s: %s terms",intFlowDirectCorrectionTermsForNUAName.Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s terms)",sinCosFlag[sc].Data()),10,0,10,"s");
1782    fNestedLoopsList->Add(fIntFlowDirectCorrectionTermsForNUA[sc]);
1783   } // end of for(Int_t sc=0;sc<2;sc++) 
1784  } // end of if(fEvaluateIntFlowNestedLoops)
1785  
1786  // nested loops for differential flow: 
1787  if(fEvaluateDiffFlowNestedLoops)
1788  {
1789   // reduced correlations:
1790   TString diffFlowDirectCorrelationsName = "fDiffFlowDirectCorrelations";
1791   diffFlowDirectCorrelationsName += fAnalysisLabel->Data();
1792   for(Int_t t=0;t<2;t++) // type: RP or POI
1793   { 
1794    for(Int_t pe=0;pe<2;pe++) // pt or eta
1795    {
1796     for(Int_t rci=0;rci<4;rci++) // reduced correlation index
1797     {
1798      // reduced correlations:
1799      fDiffFlowDirectCorrelations[t][pe][rci] = new TProfile(Form("%s, %s, %s, %s",diffFlowDirectCorrelationsName.Data(),typeFlag[t].Data(),ptEtaFlag[pe].Data(),reducedCorrelationIndex[rci].Data()),Form("%s, %s, %s, %s",diffFlowDirectCorrelationsName.Data(),typeFlag[t].Data(),ptEtaFlag[pe].Data(),reducedCorrelationIndex[rci].Data()),1,lowerPtEtaEdge[pe],upperPtEtaEdge[pe],"s");
1800      fDiffFlowDirectCorrelations[t][pe][rci]->SetXTitle(ptEtaFlag[pe].Data());
1801      fNestedLoopsList->Add(fDiffFlowDirectCorrelations[t][pe][rci]); // to be improved (add dedicated list to hold reduced correlations)
1802     } // end of for(Int_t rci=0;rci<4;rci++) // correlation index
1803    } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta 
1804   } // end of for(Int_t t=0;t<2;t++) // type: RP or POI 
1805  } // end of if(fEvaluateDiffFlowNestedLoops)
1806  // correction terms for non-uniform acceptance:
1807  TString diffFlowDirectCorrectionTermsForNUAName = "fDiffFlowDirectCorrectionTermsForNUA";
1808  diffFlowDirectCorrectionTermsForNUAName += fAnalysisLabel->Data();
1809  for(Int_t t=0;t<2;t++) // typeFlag (0 = RP, 1 = POI)
1810  { 
1811   for(Int_t pe=0;pe<2;pe++) // pt or eta
1812   {
1813    for(Int_t sc=0;sc<2;sc++) // sin or cos
1814    {
1815     for(Int_t cti=0;cti<9;cti++) // correction term index
1816     {
1817      fDiffFlowDirectCorrectionTermsForNUA[t][pe][sc][cti] = new TProfile(Form("%s, %s, %s, %s, cti = %d",diffFlowDirectCorrectionTermsForNUAName.Data(),typeFlag[t].Data(),ptEtaFlag[pe].Data(),sinCosFlag[sc].Data(),cti+1),Form("%s, %s, %s, %s, cti = %d",diffFlowDirectCorrectionTermsForNUAName.Data(),typeFlag[t].Data(),ptEtaFlag[pe].Data(),sinCosFlag[sc].Data(),cti+1),1,lowerPtEtaEdge[pe],upperPtEtaEdge[pe],"s"); 
1818      fNestedLoopsList->Add(fDiffFlowDirectCorrectionTermsForNUA[t][pe][sc][cti]);
1819     }
1820    }
1821   }
1822  } 
1823
1824 } // end of AliFlowAnalysisWithQCumulants::BookEverythingForNestedLoops()
1825
1826
1827 //================================================================================================================================
1828
1829
1830 void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
1831 {
1832  // calculate all correlations needed for integrated flow
1833  
1834  // multiplicity:
1835  Double_t dMult = (*fSMpk)(0,0);
1836  
1837  // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
1838  Double_t dReQ1n = (*fReQ)(0,0);
1839  Double_t dReQ2n = (*fReQ)(1,0);
1840  Double_t dReQ3n = (*fReQ)(2,0);
1841  Double_t dReQ4n = (*fReQ)(3,0);
1842  Double_t dImQ1n = (*fImQ)(0,0);
1843  Double_t dImQ2n = (*fImQ)(1,0);
1844  Double_t dImQ3n = (*fImQ)(2,0);
1845  Double_t dImQ4n = (*fImQ)(3,0);
1846   
1847  // real and imaginary parts of some expressions involving various combinations of Q-vectors evaluated in harmonics n, 2n, 3n and 4n:
1848  // (these expression appear in the Eqs. for the multi-particle correlations bellow)
1849  
1850  // Re[Q_{2n} Q_{n}^* Q_{n}^*]
1851  Double_t reQ2nQ1nstarQ1nstar = pow(dReQ1n,2.)*dReQ2n + 2.*dReQ1n*dImQ1n*dImQ2n - pow(dImQ1n,2.)*dReQ2n; 
1852  
1853  // Im[Q_{2n} Q_{n}^* Q_{n}^*]
1854  //Double_t imQ2nQ1nstarQ1nstar = pow(dReQ1n,2.)*dImQ2n-2.*dReQ1n*dImQ1n*dReQ2n-pow(dImQ1n,2.)*dImQ2n; 
1855  
1856  // Re[Q_{n} Q_{n} Q_{2n}^*] = Re[Q_{2n} Q_{n}^* Q_{n}^*]
1857  Double_t reQ1nQ1nQ2nstar = reQ2nQ1nstarQ1nstar; 
1858  
1859  // Re[Q_{3n} Q_{n} Q_{2n}^* Q_{2n}^*]
1860  Double_t reQ3nQ1nQ2nstarQ2nstar = (pow(dReQ2n,2.)-pow(dImQ2n,2.))*(dReQ3n*dReQ1n-dImQ3n*dImQ1n) 
1861                                  + 2.*dReQ2n*dImQ2n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);
1862
1863  // Im[Q_{3n} Q_{n} Q_{2n}^* Q_{2n}^*]                                                                  
1864  //Double_t imQ3nQ1nQ2nstarQ2nstar = calculate and implement this (deleteMe)
1865   
1866  // Re[Q_{2n} Q_{2n} Q_{3n}^* Q_{1n}^*] = Re[Q_{3n} Q_{n} Q_{2n}^* Q_{2n}^*]
1867  Double_t reQ2nQ2nQ3nstarQ1nstar = reQ3nQ1nQ2nstarQ2nstar;
1868   
1869  // Re[Q_{4n} Q_{2n}^* Q_{2n}^*]
1870  Double_t reQ4nQ2nstarQ2nstar = pow(dReQ2n,2.)*dReQ4n+2.*dReQ2n*dImQ2n*dImQ4n-pow(dImQ2n,2.)*dReQ4n;
1871
1872  // Im[Q_{4n} Q_{2n}^* Q_{2n}^*]
1873  //Double_t imQ4nQ2nstarQ2nstar = calculate and implement this (deleteMe)
1874  
1875  // Re[Q_{2n} Q_{2n} Q_{4n}^*] =  Re[Q_{4n} Q_{2n}^* Q_{2n}^*]
1876  Double_t reQ2nQ2nQ4nstar = reQ4nQ2nstarQ2nstar;
1877  
1878  // Re[Q_{4n} Q_{3n}^* Q_{n}^*]
1879  Double_t reQ4nQ3nstarQ1nstar = dReQ4n*(dReQ3n*dReQ1n-dImQ3n*dImQ1n)+dImQ4n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);
1880  
1881  // Re[Q_{3n} Q_{n} Q_{4n}^*] = Re[Q_{4n} Q_{3n}^* Q_{n}^*]
1882  Double_t reQ3nQ1nQ4nstar = reQ4nQ3nstarQ1nstar;
1883  
1884  // Im[Q_{4n} Q_{3n}^* Q_{n}^*]
1885  //Double_t imQ4nQ3nstarQ1nstar = calculate and implement this (deleteMe)
1886
1887  // Re[Q_{3n} Q_{2n}^* Q_{n}^*]
1888  Double_t reQ3nQ2nstarQ1nstar = dReQ3n*dReQ2n*dReQ1n-dReQ3n*dImQ2n*dImQ1n+dImQ3n*dReQ2n*dImQ1n
1889                               + dImQ3n*dImQ2n*dReQ1n;
1890                               
1891  // Re[Q_{2n} Q_{n} Q_{3n}^*] = Re[Q_{3n} Q_{2n}^* Q_{n}^*]
1892  Double_t reQ2nQ1nQ3nstar = reQ3nQ2nstarQ1nstar;
1893  
1894  // Im[Q_{3n} Q_{2n}^* Q_{n}^*]
1895  //Double_t imQ3nQ2nstarQ1nstar; //calculate and implement this (deleteMe)
1896  
1897  // Re[Q_{3n} Q_{n}^* Q_{n}^* Q_{n}^*]
1898  Double_t reQ3nQ1nstarQ1nstarQ1nstar = dReQ3n*pow(dReQ1n,3)-3.*dReQ1n*dReQ3n*pow(dImQ1n,2)
1899                                      + 3.*dImQ1n*dImQ3n*pow(dReQ1n,2)-dImQ3n*pow(dImQ1n,3);
1900
1901  // Im[Q_{3n} Q_{n}^* Q_{n}^* Q_{n}^*]
1902  //Double_t imQ3nQ1nstarQ1nstarQ1nstar; //calculate and implement this (deleteMe)
1903  
1904  // |Q_{2n}|^2 |Q_{n}|^2
1905  Double_t dQ2nQ1nQ2nstarQ1nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))*(pow(dReQ1n,2.)+pow(dImQ1n,2.));
1906  
1907  // Re[Q_{4n} Q_{2n}^* Q_{n}^* Q_{n}^*]
1908  Double_t reQ4nQ2nstarQ1nstarQ1nstar = (dReQ4n*dReQ2n+dImQ4n*dImQ2n)*(pow(dReQ1n,2)-pow(dImQ1n,2))
1909                                      + 2.*dReQ1n*dImQ1n*(dImQ4n*dReQ2n-dReQ4n*dImQ2n); 
1910  
1911  // Im[Q_{4n} Q_{2n}^* Q_{n}^* Q_{n}^*]
1912  //Double_t imQ4nQ2nstarQ1nstarQ1nstar; //calculate and implement this (deleteMe)
1913  
1914  // Re[Q_{2n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^*]
1915  Double_t reQ2nQ1nQ1nstarQ1nstarQ1nstar = (dReQ2n*dReQ1n-dImQ2n*dImQ1n)*(pow(dReQ1n,3)-3.*dReQ1n*pow(dImQ1n,2))
1916                                         + (dReQ2n*dImQ1n+dReQ1n*dImQ2n)*(3.*dImQ1n*pow(dReQ1n,2)-pow(dImQ1n,3));
1917
1918  // Im[Q_{2n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^*] 
1919  //Double_t imQ2nQ1nQ1nstarQ1nstarQ1nstar; //calculate and implement this (deleteMe)
1920  
1921  // Re[Q_{2n} Q_{2n} Q_{2n}^* Q_{n}^* Q_{n}^*]
1922  Double_t reQ2nQ2nQ2nstarQ1nstarQ1nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))
1923                                         * (dReQ2n*(pow(dReQ1n,2.)-pow(dImQ1n,2.)) + 2.*dImQ2n*dReQ1n*dImQ1n);
1924
1925  // Im[Q_{2n} Q_{2n} Q_{2n}^* Q_{n}^* Q_{n}^*]
1926  //Double_t imQ2nQ2nQ2nstarQ1nstarQ1nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))
1927  //                                       * (dImQ2n*(pow(dReQ1n,2.)-pow(dImQ1n,2.)) - 2.*dReQ2n*dReQ1n*dImQ1n);
1928  
1929  // Re[Q_{4n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
1930  Double_t reQ4nQ1nstarQ1nstarQ1nstarQ1nstar = pow(dReQ1n,4.)*dReQ4n-6.*pow(dReQ1n,2.)*dReQ4n*pow(dImQ1n,2.)
1931                                             + pow(dImQ1n,4.)*dReQ4n+4.*pow(dReQ1n,3.)*dImQ1n*dImQ4n
1932                                             - 4.*pow(dImQ1n,3.)*dReQ1n*dImQ4n;
1933                                             
1934  // Im[Q_{4n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
1935  //Double_t imQ4nQ1nstarQ1nstarQ1nstarQ1nstar = pow(dReQ1n,4.)*dImQ4n-6.*pow(dReQ1n,2.)*dImQ4n*pow(dImQ1n,2.)
1936  //                                           + pow(dImQ1n,4.)*dImQ4n+4.*pow(dImQ1n,3.)*dReQ1n*dReQ4n
1937  //                                           - 4.*pow(dReQ1n,3.)*dImQ1n*dReQ4n;
1938  
1939  // Re[Q_{3n} Q_{n} Q_{2n}^* Q_{n}^* Q_{n}^*]
1940  Double_t reQ3nQ1nQ2nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)+pow(dImQ1n,2.))
1941                                         * (dReQ1n*dReQ2n*dReQ3n-dReQ3n*dImQ1n*dImQ2n+dReQ2n*dImQ1n*dImQ3n+dReQ1n*dImQ2n*dImQ3n);
1942  
1943  // Im[Q_{3n} Q_{n} Q_{2n}^* Q_{n}^* Q_{n}^*]
1944  //Double_t imQ3nQ1nQ2nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)+pow(dImQ1n,2.))
1945  //                                       * (-dReQ2n*dReQ3n*dImQ1n-dReQ1n*dReQ3n*dImQ2n+dReQ1n*dReQ2n*dImQ3n-dImQ1n*dImQ2n*dImQ3n);
1946  
1947  
1948  // Re[Q_{2n} Q_{2n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
1949  Double_t reQ2nQ2nQ1nstarQ1nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)*dReQ2n-2.*dReQ1n*dReQ2n*dImQ1n-dReQ2n*pow(dImQ1n,2.)
1950                                                + dImQ2n*pow(dReQ1n,2.)+2.*dReQ1n*dImQ1n*dImQ2n-pow(dImQ1n,2.)*dImQ2n)
1951                                                * (pow(dReQ1n,2.)*dReQ2n+2.*dReQ1n*dReQ2n*dImQ1n-dReQ2n*pow(dImQ1n,2.)
1952                                                - dImQ2n*pow(dReQ1n,2.)+2.*dReQ1n*dImQ1n*dImQ2n+pow(dImQ1n,2.)*dImQ2n);
1953  
1954  // Im[Q_{2n} Q_{2n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
1955  //Double_t imQ2nQ2nQ1nstarQ1nstarQ1nstarQ1nstar = 2.*(pow(dReQ1n,2.)*dReQ2n-dReQ2n*pow(dImQ1n,2.)
1956  //                                              + 2.*dReQ1n*dImQ1n*dImQ2n)*(pow(dReQ1n,2.)*dImQ2n
1957  //                                              - 2.*dReQ1n*dImQ1n*dReQ2n-pow(dImQ1n,2.)*dImQ2n);
1958  
1959  // Re[Q_{3n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
1960  Double_t reQ3nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)+pow(dImQ1n,2.))
1961                                                * (pow(dReQ1n,3.)*dReQ3n-3.*dReQ1n*dReQ3n*pow(dImQ1n,2.)
1962                                                + 3.*pow(dReQ1n,2.)*dImQ1n*dImQ3n-pow(dImQ1n,3.)*dImQ3n);
1963   
1964  // Im[Q_{3n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]                                                                                           
1965  //Double_t imQ3nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)+pow(dImQ1n,2.))
1966  //                                              * (pow(dImQ1n,3.)*dReQ3n-3.*dImQ1n*dReQ3n*pow(dReQ1n,2.)
1967  //                                              - 3.*pow(dImQ1n,2.)*dReQ1n*dImQ3n+pow(dReQ1n,3.)*dImQ3n);
1968  
1969  // |Q_{2n}|^2 |Q_{n}|^4
1970  Double_t dQ2nQ1nQ1nQ2nstarQ1nstarQ1nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))*pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.);
1971  
1972  // Re[Q_{2n} Q_{n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]
1973  Double_t reQ2nQ1nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)
1974                                                   * (pow(dReQ1n,2.)*dReQ2n-dReQ2n*pow(dImQ1n,2.)
1975                                                   + 2.*dReQ1n*dImQ1n*dImQ2n);
1976                                                   
1977  // Im[Q_{2n} Q_{n} Q_{n} Q_{n}^* Q_{n}^* Q_{n}^* Q_{n}^*]                                                  
1978  //Double_t imQ2nQ1nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)
1979  //                                                 * (pow(dReQ1n,2.)*dImQ2n-dImQ2n*pow(dImQ1n,2.)
1980  //                                                 - 2.*dReQ1n*dReQ2n*dImQ1n);
1981  
1982   
1983  
1984        
1985  //                                        **************************************
1986  //                                        **** multi-particle correlations: ****
1987  //                                        **************************************
1988  //
1989  // Remark 1: multi-particle correlations calculated with non-weighted Q-vectors are stored in 1D profile fQCorrelations[0]. // to be improved (wrong profiles)
1990  // Remark 2: binning of fQCorrelations[0] is organized as follows: // to be improved (wrong profiles)
1991  // --------------------------------------------------------------------------------------------------------------------
1992  //  1st bin: <2>_{1n|1n} = two1n1n = cos(n*(phi1-phi2))>
1993  //  2nd bin: <2>_{2n|2n} = two2n2n = cos(2n*(phi1-phi2))>
1994  //  3rd bin: <2>_{3n|3n} = two3n3n = cos(3n*(phi1-phi2))> 
1995  //  4th bin: <2>_{4n|4n} = two4n4n = cos(4n*(phi1-phi2))>
1996  //  5th bin:           ----  EMPTY ----
1997  //  6th bin: <3>_{2n|1n,1n} = three2n1n1n = <cos(n*(2.*phi1-phi2-phi3))>
1998  //  7th bin: <3>_{3n|2n,1n} = three3n2n1n = <cos(n*(3.*phi1-2.*phi2-phi3))>
1999  //  8th bin: <3>_{4n|2n,2n} = three4n2n2n = <cos(n*(4.*phi1-2.*phi2-2.*phi3))>
2000  //  9th bin: <3>_{4n|3n,1n} = three4n3n1n = <cos(n*(4.*phi1-3.*phi2-phi3))>
2001  // 10th bin:           ----  EMPTY ----
2002  // 11th bin: <4>_{1n,1n|1n,1n} = four1n1n1n1n = <cos(n*(phi1+phi2-phi3-phi4))>
2003  // 12th bin: <4>_{2n,1n|2n,1n} = four2n1n2n1n = <cos(2.*n*(phi1+phi2-phi3-phi4))>
2004  // 13th bin: <4>_{2n,2n|2n,2n} = four2n2n2n2n = <cos(n*(2.*phi1+phi2-2.*phi3-phi4))>
2005  // 14th bin: <4>_{3n|1n,1n,1n} = four3n1n1n1n = <cos(n*(3.*phi1-phi2-phi3-phi4))> 
2006  // 15th bin: <4>_{3n,1n|3n,1n} = four3n1n3n1n = <cos(n*(4.*phi1-2.*phi2-phi3-phi4))>
2007  // 16th bin: <4>_{3n,1n|2n,2n} = four3n1n2n2n = <cos(n*(3.*phi1+phi2-2.*phi3-2.*phi4))>
2008  // 17th bin: <4>_{4n|2n,1n,1n} = four4n2n1n1n = <cos(n*(3.*phi1+phi2-3.*phi3-phi4))> 
2009  // 18th bin:           ----  EMPTY ----
2010  // 19th bin: <5>_{2n|1n,1n,1n,1n} = five2n1n1n1n1n = <cos(n*(2.*phi1+phi2-phi3-phi4-phi5))>
2011  // 20th bin: <5>_{2n,2n|2n,1n,1n} = five2n2n2n1n1n = <cos(n*(2.*phi1+2.*phi2-2.*phi3-phi4-phi5))>
2012  // 21st bin: <5>_{3n,1n|2n,1n,1n} = five3n1n2n1n1n = <cos(n*(3.*phi1+phi2-2.*phi3-phi4-phi5))>
2013  // 22nd bin: <5>_{4n|1n,1n,1n,1n} = five4n1n1n1n1n = <cos(n*(4.*phi1-phi2-phi3-phi4-phi5))>
2014  // 23rd bin:           ----  EMPTY ----
2015  // 24th bin: <6>_{1n,1n,1n|1n,1n,1n} = six1n1n1n1n1n1n = <cos(n*(phi1+phi2+phi3-phi4-phi5-phi6))>
2016  // 25th bin: <6>_{2n,1n,1n|2n,1n,1n} = six2n1n1n2n1n1n = <cos(n*(2.*phi1+2.*phi2-phi3-phi4-phi5-phi6))>
2017  // 26th bin: <6>_{2n,2n|1n,1n,1n,1n} = six2n2n1n1n1n1n = <cos(n*(3.*phi1+phi2-phi3-phi4-phi5-phi6))>
2018  // 27th bin: <6>_{3n,1n|1n,1n,1n,1n} = six3n1n1n1n1n1n = <cos(n*(2.*phi1+phi2+phi3-2.*phi4-phi5-phi6))>
2019  // 28th bin:           ----  EMPTY ----
2020  // 29th bin: <7>_{2n,1n,1n|1n,1n,1n,1n} = seven2n1n1n1n1n1n1n =  <cos(n*(2.*phi1+phi2+phi3-phi4-phi5-phi6-phi7))>
2021  // 30th bin:           ----  EMPTY ----
2022  // 31st bin: <8>_{1n,1n,1n,1n|1n,1n,1n,1n} = eight1n1n1n1n1n1n1n1n = <cos(n*(phi1+phi2+phi3+phi4-phi5-phi6-phi7-phi8))>
2023  // --------------------------------------------------------------------------------------------------------------------
2024     
2025  // 2-particle:
2026  Double_t two1n1n = 0.; // <cos(n*(phi1-phi2))>
2027  Double_t two2n2n = 0.; // <cos(2n*(phi1-phi2))>
2028  Double_t two3n3n = 0.; // <cos(3n*(phi1-phi2))>
2029  Double_t two4n4n = 0.; // <cos(4n*(phi1-phi2))>
2030  
2031  if(dMult>1)
2032  {
2033   two1n1n = (pow(dReQ1n,2.)+pow(dImQ1n,2.)-dMult)/(dMult*(dMult-1.)); 
2034   two2n2n = (pow(dReQ2n,2.)+pow(dImQ2n,2.)-dMult)/(dMult*(dMult-1.)); 
2035   two3n3n = (pow(dReQ3n,2.)+pow(dImQ3n,2.)-dMult)/(dMult*(dMult-1.)); 
2036   two4n4n = (pow(dReQ4n,2.)+pow(dImQ4n,2.)-dMult)/(dMult*(dMult-1.)); 
2037   
2038   // average 2-particle correlations for single event: 
2039   fIntFlowCorrelationsAllEBE->SetBinContent(1,two1n1n);
2040   fIntFlowCorrelationsAllEBE->SetBinContent(2,two2n2n);
2041   fIntFlowCorrelationsAllEBE->SetBinContent(3,two3n3n);
2042   fIntFlowCorrelationsAllEBE->SetBinContent(4,two4n4n);
2043           
2044   // average 2-particle correlations for all events:      
2045   fIntFlowCorrelationsAllPro->Fill(0.5,two1n1n,dMult*(dMult-1.));  
2046   fIntFlowCorrelationsAllPro->Fill(1.5,two2n2n,dMult*(dMult-1.)); 
2047   fIntFlowCorrelationsAllPro->Fill(2.5,two3n3n,dMult*(dMult-1.)); 
2048   fIntFlowCorrelationsAllPro->Fill(3.5,two4n4n,dMult*(dMult-1.)); 
2049   
2050   // store separetately <2> (to be improved: do I really need this?)
2051   fIntFlowCorrelationsEBE->SetBinContent(1,two1n1n); // <2>
2052   fIntFlowEventWeightsForCorrelationsEBE->SetBinContent(1,dMult*(dMult-1.)); // eW_<2>
2053   fIntFlowCorrelationsPro->Fill(0.5,two1n1n,dMult*(dMult-1.));
2054   
2055   // distribution of <cos(n*(phi1-phi2))>:
2056   //f2pDistribution->Fill(two1n1n,dMult*(dMult-1.)); 
2057  } // end of if(dMult>1)
2058  
2059  // 3-particle:
2060  Double_t three2n1n1n = 0.; // <cos(n*(2.*phi1-phi2-phi3))>
2061  Double_t three3n2n1n = 0.; // <cos(n*(3.*phi1-2.*phi2-phi3))>
2062  Double_t three4n2n2n = 0.; // <cos(n*(4.*phi1-2.*phi2-2.*phi3))>
2063  Double_t three4n3n1n = 0.; // <cos(n*(4.*phi1-3.*phi2-phi3))>
2064  
2065  if(dMult>2)
2066  {
2067   three2n1n1n = (reQ2nQ1nstarQ1nstar-2.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
2068               - (pow(dReQ2n,2.)+pow(dImQ2n,2.))+2.*dMult)
2069               / (dMult*(dMult-1.)*(dMult-2.));              
2070   three3n2n1n = (reQ3nQ2nstarQ1nstar-(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2071               - (pow(dReQ2n,2.)+pow(dImQ2n,2.))
2072               - (pow(dReQ1n,2.)+pow(dImQ1n,2.))+2.*dMult)
2073               / (dMult*(dMult-1.)*(dMult-2.));
2074   three4n2n2n = (reQ4nQ2nstarQ2nstar-2.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
2075               - (pow(dReQ4n,2.)+pow(dImQ4n,2.))+2.*dMult)
2076               / (dMult*(dMult-1.)*(dMult-2.)); 
2077   three4n3n1n = (reQ4nQ3nstarQ1nstar-(pow(dReQ4n,2.)+pow(dImQ4n,2.))
2078               - (pow(dReQ3n,2.)+pow(dImQ3n,2.))
2079               - (pow(dReQ1n,2.)+pow(dImQ1n,2.))+2.*dMult)
2080               / (dMult*(dMult-1.)*(dMult-2.)); 
2081               
2082   // average 3-particle correlations for single event: 
2083   fIntFlowCorrelationsAllEBE->SetBinContent(6,three2n1n1n);
2084   fIntFlowCorrelationsAllEBE->SetBinContent(7,three3n2n1n);
2085   fIntFlowCorrelationsAllEBE->SetBinContent(8,three4n2n2n);
2086   fIntFlowCorrelationsAllEBE->SetBinContent(9,three4n3n1n);
2087         
2088   // average 3-particle correlations for all events:                
2089   fIntFlowCorrelationsAllPro->Fill(5.5,three2n1n1n,dMult*(dMult-1.)*(dMult-2.)); 
2090   fIntFlowCorrelationsAllPro->Fill(6.5,three3n2n1n,dMult*(dMult-1.)*(dMult-2.));
2091   fIntFlowCorrelationsAllPro->Fill(7.5,three4n2n2n,dMult*(dMult-1.)*(dMult-2.)); 
2092   fIntFlowCorrelationsAllPro->Fill(8.5,three4n3n1n,dMult*(dMult-1.)*(dMult-2.));    
2093  } // end of if(dMult>2)
2094  
2095  // 4-particle:
2096  Double_t four1n1n1n1n = 0.; // <cos(n*(phi1+phi2-phi3-phi4))>
2097  Double_t four2n2n2n2n = 0.; // <cos(2.*n*(phi1+phi2-phi3-phi4))>
2098  Double_t four2n1n2n1n = 0.; // <cos(n*(2.*phi1+phi2-2.*phi3-phi4))> 
2099  Double_t four3n1n1n1n = 0.; // <cos(n*(3.*phi1-phi2-phi3-phi4))> 
2100  Double_t four4n2n1n1n = 0.; // <cos(n*(4.*phi1-2.*phi2-phi3-phi4))> 
2101  Double_t four3n1n2n2n = 0.; // <cos(n*(3.*phi1+phi2-2.*phi3-2.*phi4))> 
2102  Double_t four3n1n3n1n = 0.; // <cos(n*(3.*phi1+phi2-3.*phi3-phi4))>   
2103  
2104  if(dMult>3)
2105  {
2106   four1n1n1n1n = (2.*dMult*(dMult-3.)+pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)-4.*(dMult-2.)*(pow(dReQ1n,2.)
2107                + pow(dImQ1n,2.))-2.*reQ2nQ1nstarQ1nstar+(pow(dReQ2n,2.)+pow(dImQ2n,2.)))
2108                / (dMult*(dMult-1)*(dMult-2.)*(dMult-3.));     
2109   four2n2n2n2n = (2.*dMult*(dMult-3.)+pow((pow(dReQ2n,2.)+pow(dImQ2n,2.)),2.)-4.*(dMult-2.)*(pow(dReQ2n,2.)
2110                + pow(dImQ2n,2.))-2.*reQ4nQ2nstarQ2nstar+(pow(dReQ4n,2.)+pow(dImQ4n,2.)))
2111                / (dMult*(dMult-1)*(dMult-2.)*(dMult-3.));
2112   four2n1n2n1n = (dQ2nQ1nQ2nstarQ1nstar-2.*reQ3nQ2nstarQ1nstar-2.*reQ2nQ1nstarQ1nstar)
2113                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2114                - ((dMult-5.)*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
2115                + (dMult-4.)*(pow(dReQ2n,2.)+pow(dImQ2n,2.))-(pow(dReQ3n,2.)+pow(dImQ3n,2.)))
2116                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2117                + (dMult-6.)/((dMult-1.)*(dMult-2.)*(dMult-3.));
2118   four3n1n1n1n = (reQ3nQ1nstarQ1nstarQ1nstar-3.*reQ3nQ2nstarQ1nstar-3.*reQ2nQ1nstarQ1nstar)
2119                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2120                + (2.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))+3.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
2121                + 6.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))-6.*dMult)
2122                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2123   four4n2n1n1n = (reQ4nQ2nstarQ1nstarQ1nstar-2.*reQ4nQ3nstarQ1nstar-reQ4nQ2nstarQ2nstar-2.*reQ3nQ2nstarQ1nstar)
2124                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2125                - (reQ2nQ1nstarQ1nstar-2.*(pow(dReQ4n,2.)+pow(dImQ4n,2.))-2.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2126                - 3.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))-4.*(pow(dReQ1n,2.)+pow(dImQ1n,2.)))
2127                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2128                - 6./((dMult-1.)*(dMult-2.)*(dMult-3.));
2129   four3n1n2n2n = (reQ3nQ1nQ2nstarQ2nstar-reQ4nQ2nstarQ2nstar-reQ3nQ1nQ4nstar-2.*reQ3nQ2nstarQ1nstar)
2130                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2131                - (2.*reQ1nQ1nQ2nstar-(pow(dReQ4n,2.)+pow(dImQ4n,2.))-2.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2132                - 4.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))-4.*(pow(dReQ1n,2.)+pow(dImQ1n,2.)))
2133                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2134                - 6./((dMult-1.)*(dMult-2.)*(dMult-3.)); 
2135   four3n1n3n1n = ((pow(dReQ3n,2.)+pow(dImQ3n,2.))*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
2136                - 2.*reQ4nQ3nstarQ1nstar-2.*reQ3nQ2nstarQ1nstar)
2137                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2138                + ((pow(dReQ4n,2.)+pow(dImQ4n,2.))-(dMult-4.)*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2139                + (pow(dReQ2n,2.)+pow(dImQ2n,2.))-(dMult-4.)*(pow(dReQ1n,2.)+pow(dImQ1n,2.)))
2140                / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.))
2141                + (dMult-6.)/((dMult-1.)*(dMult-2.)*(dMult-3.));
2142                
2143   // average 4-particle correlations for single event: 
2144   fIntFlowCorrelationsAllEBE->SetBinContent(11,four1n1n1n1n);
2145   fIntFlowCorrelationsAllEBE->SetBinContent(12,four2n1n2n1n);
2146   fIntFlowCorrelationsAllEBE->SetBinContent(13,four2n2n2n2n);
2147   fIntFlowCorrelationsAllEBE->SetBinContent(14,four3n1n1n1n);
2148   fIntFlowCorrelationsAllEBE->SetBinContent(15,four3n1n3n1n);
2149   fIntFlowCorrelationsAllEBE->SetBinContent(16,four3n1n2n2n);
2150   fIntFlowCorrelationsAllEBE->SetBinContent(17,four4n2n1n1n);
2151         
2152   // average 4-particle correlations for all events:                
2153   fIntFlowCorrelationsAllPro->Fill(10.5,four1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2154   fIntFlowCorrelationsAllPro->Fill(11.5,four2n1n2n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2155   fIntFlowCorrelationsAllPro->Fill(12.5,four2n2n2n2n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2156   fIntFlowCorrelationsAllPro->Fill(13.5,four3n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2157   fIntFlowCorrelationsAllPro->Fill(14.5,four3n1n3n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2158   fIntFlowCorrelationsAllPro->Fill(15.5,four3n1n2n2n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));  
2159   fIntFlowCorrelationsAllPro->Fill(16.5,four4n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)); 
2160   
2161   // store separetately <4> (to be improved: do I really need this?)
2162   fIntFlowCorrelationsEBE->SetBinContent(2,four1n1n1n1n); // <4>
2163   fIntFlowEventWeightsForCorrelationsEBE->SetBinContent(2,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)); // eW_<4>
2164   fIntFlowCorrelationsPro->Fill(1.5,four1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2165   
2166   // distribution of <cos(n*(phi1+phi2-phi3-phi4))>
2167   //f4pDistribution->Fill(four1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.));
2168   
2169  } // end of if(dMult>3)
2170
2171  // 5-particle:
2172  Double_t five2n1n1n1n1n = 0.; // <cos(n*(2.*phi1+phi2-phi3-phi4-phi5))>
2173  Double_t five2n2n2n1n1n = 0.; // <cos(n*(2.*phi1+2.*phi2-2.*phi3-phi4-phi5))>
2174  Double_t five3n1n2n1n1n = 0.; // <cos(n*(3.*phi1+phi2-2.*phi3-phi4-phi5))>
2175  Double_t five4n1n1n1n1n = 0.; // <cos(n*(4.*phi1-phi2-phi3-phi4-phi5))>
2176  
2177  if(dMult>4)
2178  {
2179   five2n1n1n1n1n = (reQ2nQ1nQ1nstarQ1nstarQ1nstar-reQ3nQ1nstarQ1nstarQ1nstar+6.*reQ3nQ2nstarQ1nstar)
2180                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2181                  - (reQ2nQ1nQ3nstar+3.*(dMult-6.)*reQ2nQ1nstarQ1nstar+3.*reQ1nQ1nQ2nstar)
2182                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2183                  - (2.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2184                  + 3.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))*(pow(dReQ2n,2.)+pow(dImQ2n,2.))     
2185                  - 3.*(dMult-4.)*(pow(dReQ2n,2.)+pow(dImQ2n,2.)))
2186                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2187                  - 3.*(pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)
2188                  - 2.*(2*dMult-5.)*(pow(dReQ1n,2.)+pow(dImQ1n,2.))+2.*dMult*(dMult-4.))
2189                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2190                  
2191   five2n2n2n1n1n = (reQ2nQ2nQ2nstarQ1nstarQ1nstar-reQ4nQ2nstarQ1nstarQ1nstar-2.*reQ2nQ2nQ3nstarQ1nstar)
2192                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2193                  + 2.*(reQ4nQ2nstarQ2nstar+4.*reQ3nQ2nstarQ1nstar+reQ3nQ1nQ4nstar)
2194                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2195                  + (reQ2nQ2nQ4nstar-2.*(dMult-5.)*reQ2nQ1nstarQ1nstar+2.*reQ1nQ1nQ2nstar)
2196                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2197                  - (2.*(pow(dReQ4n,2.)+pow(dImQ4n,2.))+4.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2198                  + 1.*pow((pow(dReQ2n,2.)+pow(dImQ2n,2.)),2.)
2199                  - 2.*(3.*dMult-10.)*(pow(dReQ2n,2.)+pow(dImQ2n,2.)))
2200                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2201                  - (4.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
2202                  - 4.*(dMult-5.)*(pow(dReQ1n,2.)+pow(dImQ1n,2.))+4.*dMult*(dMult-6.))
2203                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)); 
2204
2205   five4n1n1n1n1n = (reQ4nQ1nstarQ1nstarQ1nstarQ1nstar-6.*reQ4nQ2nstarQ1nstarQ1nstar-4.*reQ3nQ1nstarQ1nstarQ1nstar)
2206                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2207                  + (8.*reQ4nQ3nstarQ1nstar+3.*reQ4nQ2nstarQ2nstar+12.*reQ3nQ2nstarQ1nstar+12.*reQ2nQ1nstarQ1nstar)
2208                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2209                  - (6.*(pow(dReQ4n,2.)+pow(dImQ4n,2.))+8.*(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2210                  + 12.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))+24.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))-24.*dMult)
2211                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2212   
2213   five3n1n2n1n1n = (reQ3nQ1nQ2nstarQ1nstarQ1nstar-reQ4nQ2nstarQ1nstarQ1nstar-reQ3nQ1nstarQ1nstarQ1nstar)
2214                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2215                  - (reQ3nQ1nQ2nstarQ2nstar-3.*reQ4nQ3nstarQ1nstar-reQ4nQ2nstarQ2nstar)
2216                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2217                  - ((2.*dMult-13.)*reQ3nQ2nstarQ1nstar-reQ3nQ1nQ4nstar-9.*reQ2nQ1nstarQ1nstar)
2218                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2219                  - (2.*reQ1nQ1nQ2nstar+2.*(pow(dReQ4n,2.)+pow(dImQ4n,2.))
2220                  - 2.*(dMult-5.)*(pow(dReQ3n,2.)+pow(dImQ3n,2.))+2.*(pow(dReQ3n,2.)
2221                  + pow(dImQ3n,2.))*(pow(dReQ1n,2.)+pow(dImQ1n,2.)))
2222                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2223                  + (2.*(dMult-6.)*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
2224                  - 2.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
2225                  - pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)
2226                  + 2.*(3.*dMult-11.)*(pow(dReQ1n,2.)+pow(dImQ1n,2.)))
2227                  / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.))
2228                  - 4.*(dMult-6.)/((dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2229                  
2230   // average 5-particle correlations for single event: 
2231   fIntFlowCorrelationsAllEBE->SetBinContent(19,five2n1n1n1n1n);
2232   fIntFlowCorrelationsAllEBE->SetBinContent(20,five2n2n2n1n1n);
2233   fIntFlowCorrelationsAllEBE->SetBinContent(21,five3n1n2n1n1n);
2234   fIntFlowCorrelationsAllEBE->SetBinContent(22,five4n1n1n1n1n);
2235         
2236   // average 5-particle correlations for all events:                         
2237   fIntFlowCorrelationsAllPro->Fill(18.5,five2n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)); 
2238   fIntFlowCorrelationsAllPro->Fill(19.5,five2n2n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2239   fIntFlowCorrelationsAllPro->Fill(20.5,five3n1n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2240   fIntFlowCorrelationsAllPro->Fill(21.5,five4n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.));
2241  } // end of if(dMult>4)
2242     
2243  // 6-particle:
2244  Double_t six1n1n1n1n1n1n = 0.; // <cos(n*(phi1+phi2+phi3-phi4-phi5-phi6))>
2245  Double_t six2n2n1n1n1n1n = 0.; // <cos(n*(2.*phi1+2.*phi2-phi3-phi4-phi5-phi6))>
2246  Double_t six3n1n1n1n1n1n = 0.; // <cos(n*(3.*phi1+phi2-phi3-phi4-phi5-phi6))>
2247  Double_t six2n1n1n2n1n1n = 0.; // <cos(n*(2.*phi1+phi2+phi3-2.*phi4-phi5-phi6))>
2248  
2249  if(dMult>5)
2250  {
2251   six1n1n1n1n1n1n = (pow(pow(dReQ1n,2.)+pow(dImQ1n,2.),3.)+9.*dQ2nQ1nQ2nstarQ1nstar-6.*reQ2nQ1nQ1nstarQ1nstarQ1nstar)
2252                   / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.))
2253                   + 4.*(reQ3nQ1nstarQ1nstarQ1nstar-3.*reQ3nQ2nstarQ1nstar)
2254                   / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.))
2255                   + 2.*(9.*(dMult-4.)*reQ2nQ1nstarQ1nstar+2.*(pow(dReQ3n,2.)+pow(dImQ3n,2.)))
2256                   / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.))
2257                   - 9.*(pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)+(pow(dReQ2n,2.)+pow(dImQ2n,2.)))
2258                   / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-5.))
2259                   + (18.*(pow(dReQ1n,2.)+pow(dImQ1n,2.)))
2260                   / (dMult*(dMult-1)*(dMult-3)*(dMult-4))
2261                   - 6./((dMult-1.)*(dMult-2.)*(dMult-3.));
2262                   
2263   six2n1n1n2n1n1n = (dQ2nQ1nQ1nQ2nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)
2264                   * (2.*five2n2n2n1n1n+4.*five2n1n1n1n1n+4.*five3n1n2n1n1n+4.*four2n1n2n1n+1.*four1n1n1n1n)
2265                   - dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(4.*four1n1n1n1n+4.*two1n1n
2266                   + 2.*three2n1n1n+2.*three2n1n1n+4.*four3n1n1n1n+8.*three2n1n1n+2.*four4n2n1n1n
2267                   + 4.*four2n1n2n1n+2.*two2n2n+8.*four2n1n2n1n+4.*four3n1n3n1n+8.*three3n2n1n
2268                   + 4.*four3n1n2n2n+4.*four1n1n1n1n+4.*four2n1n2n1n+1.*four2n2n2n2n)
2269                   - dMult*(dMult-1.)*(dMult-2.)*(2.*three2n1n1n+8.*two1n1n+4.*two1n1n+2.
2270                   + 4.*two1n1n+4.*three2n1n1n+2.*two2n2n+4.*three2n1n1n+8.*three3n2n1n
2271                   + 8.*two2n2n+4.*three4n3n1n+4.*two3n3n+4.*three3n2n1n+4.*two1n1n
2272                   + 8.*three2n1n1n+4.*two1n1n+4.*three3n2n1n+4.*three2n1n1n+2.*two2n2n
2273                   + 4.*three3n2n1n+2.*three4n2n2n)-dMult*(dMult-1.)
2274                   * (4.*two1n1n+4.+4.*two1n1n+2.*two2n2n+1.+4.*two1n1n+4.*two2n2n+4.*two3n3n
2275                   + 1.+2.*two2n2n+1.*two4n4n)-dMult)
2276                   / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); // to be improved (direct formula needed)
2277  
2278   six2n2n1n1n1n1n = (reQ2nQ2nQ1nstarQ1nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)
2279                   * (five4n1n1n1n1n+8.*five2n1n1n1n1n+6.*five2n2n2n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)
2280                   * (4.*four3n1n1n1n+6.*four4n2n1n1n+12.*three2n1n1n+12.*four1n1n1n1n+24.*four2n1n2n1n
2281                   + 4.*four3n1n2n2n+3.*four2n2n2n2n)-dMult*(dMult-1.)*(dMult-2.)*(6.*three2n1n1n+12.*three3n2n1n
2282                   + 4.*three4n3n1n+3.*three4n2n2n+8.*three2n1n1n+24.*two1n1n+12.*two2n2n+12.*three2n1n1n+8.*three3n2n1n
2283                   + 1.*three4n2n2n)-dMult*(dMult-1.)*(4.*two1n1n+6.*two2n2n+4.*two3n3n+1.*two4n4n+2.*two2n2n+8.*two1n1n+6.)-dMult)
2284                   / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); // to be improved (direct formula needed)
2285    
2286   six3n1n1n1n1n1n = (reQ3nQ1nQ1nstarQ1nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)
2287                   * (five4n1n1n1n1n+4.*five2n1n1n1n1n+6.*five3n1n2n1n1n+4.*four3n1n1n1n)
2288                   - dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(4.*four3n1n1n1n+6.*four4n2n1n1n+6.*four1n1n1n1n
2289                   + 12.*three2n1n1n+12.*four2n1n2n1n+6.*four3n1n1n1n+12.*three3n2n1n+4.*four3n1n3n1n+3.*four3n1n2n2n)
2290                   - dMult*(dMult-1.)*(dMult-2.)*(6.*three2n1n1n+12.*three3n2n1n+4.*three4n3n1n+3.*three4n2n2n+4.*two1n1n
2291                   + 12.*two1n1n+6.*three2n1n1n+12.*three2n1n1n+4.*three3n2n1n+12.*two2n2n+4.*three3n2n1n+4.*two3n3n+1.*three4n3n1n
2292                   + 6.*three3n2n1n)-dMult*(dMult-1.)*(4.*two1n1n+6.*two2n2n+4.*two3n3n+1.*two4n4n+1.*two1n1n+4.+6.*two1n1n+4.*two2n2n
2293                   + 1.*two3n3n)-dMult)/(dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); // to be improved (direct formula needed)
2294                                  
2295   // average 6-particle correlations for single event: 
2296   fIntFlowCorrelationsAllEBE->SetBinContent(24,six1n1n1n1n1n1n);
2297   fIntFlowCorrelationsAllEBE->SetBinContent(25,six2n1n1n2n1n1n);
2298   fIntFlowCorrelationsAllEBE->SetBinContent(26,six2n2n1n1n1n1n);
2299   fIntFlowCorrelationsAllEBE->SetBinContent(27,six3n1n1n1n1n1n);
2300         
2301   // average 6-particle correlations for all events:         
2302   fIntFlowCorrelationsAllPro->Fill(23.5,six1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
2303   fIntFlowCorrelationsAllPro->Fill(24.5,six2n1n1n2n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
2304   fIntFlowCorrelationsAllPro->Fill(25.5,six2n2n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));
2305   fIntFlowCorrelationsAllPro->Fill(26.5,six3n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
2306
2307   // store separetately <6> (to be improved: do I really need this?)
2308   fIntFlowCorrelationsEBE->SetBinContent(3,six1n1n1n1n1n1n); // <6>
2309   fIntFlowEventWeightsForCorrelationsEBE->SetBinContent(3,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); // eW_<6>
2310   fIntFlowCorrelationsPro->Fill(2.5,six1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.));
2311  
2312   // distribution of <cos(n*(phi1+phi2+phi3-phi4-phi5-phi6))>
2313   //f6pDistribution->Fill(six1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)); 
2314  } // end of if(dMult>5)
2315  
2316  // 7-particle:
2317  Double_t seven2n1n1n1n1n1n1n = 0.; // <cos(n*(2.*phi1+phi2+phi3-phi4-phi5-phi6-phi7))>
2318  
2319  if(dMult>6)
2320  {
2321   seven2n1n1n1n1n1n1n = (reQ2nQ1nQ1nQ1nstarQ1nstarQ1nstarQ1nstar-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)
2322                       * (2.*six3n1n1n1n1n1n+4.*six1n1n1n1n1n1n+1.*six2n2n1n1n1n1n+6.*six2n1n1n2n1n1n+8.*five2n1n1n1n1n)
2323                       - dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(1.*five4n1n1n1n1n +8.*five2n1n1n1n1n+8.*four3n1n1n1n
2324                       + 12.*five3n1n2n1n1n+4.*five2n1n1n1n1n+3.*five2n2n2n1n1n+6.*five2n2n2n1n1n+6.*four1n1n1n1n+24.*four1n1n1n1n
2325                       + 12.*five2n1n1n1n1n+12.*five2n1n1n1n1n+12.*three2n1n1n+24.*four2n1n2n1n+4.*five3n1n2n1n1n+4.*five2n1n1n1n1n)
2326                       - dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(4.*four3n1n1n1n+6.*four4n2n1n1n+12.*four1n1n1n1n+24.*three2n1n1n
2327                       + 24.*four2n1n2n1n+12.*four3n1n1n1n+24.*three3n2n1n+8.*four3n1n3n1n+6.*four3n1n2n2n+6.*three2n1n1n+12.*four1n1n1n1n
2328                       + 12.*four2n1n2n1n+6.*three2n1n1n+12.*four2n1n2n1n+4.*four3n1n2n2n+3.*four2n2n2n2n+4.*four1n1n1n1n+6.*three2n1n1n
2329                       + 24.*two1n1n+24.*four1n1n1n1n+4.*four3n1n1n1n+24.*two1n1n+24.*three2n1n1n+12.*two2n2n+24.*three2n1n1n+12.*four2n1n2n1n
2330                       + 8.*three3n2n1n+8.*four2n1n2n1n+1.*four4n2n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(6.*three2n1n1n+1.*three2n1n1n+8.*two1n1n
2331                       + 12.*three3n2n1n+24.*two1n1n+12.*three2n1n1n+4.*three2n1n1n+8.*two1n1n+4.*three4n3n1n+24.*three2n1n1n+8.*three3n2n1n
2332                       + 12.*two1n1n+12.*two1n1n+3.*three4n2n2n+24.*two2n2n+6.*two2n2n+12.+12.*three3n2n1n+8.*two3n3n+12.*three2n1n1n+24.*two1n1n
2333                       + 4.*three3n2n1n+8.*three3n2n1n+2.*three4n3n1n+12.*two1n1n+8.*three2n1n1n+4.*three2n1n1n+2.*three3n2n1n+6.*two2n2n+8.*two2n2n
2334                       + 1.*three4n2n2n+4.*three3n2n1n+6.*three2n1n1n)-dMult*(dMult-1.)*(4.*two1n1n+2.*two1n1n+6.*two2n2n+8.+1.*two2n2n+4.*two3n3n
2335                       + 12.*two1n1n+4.*two1n1n+1.*two4n4n+8.*two2n2n+6.+2.*two3n3n+4.*two1n1n+1.*two2n2n)-dMult)
2336                       / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)); // to be improved (direct formula needed)
2337         
2338   // average 7-particle correlations for single event: 
2339   fIntFlowCorrelationsAllEBE->SetBinContent(29,seven2n1n1n1n1n1n1n);
2340        
2341   // average 7-particle correlations for all events:                      
2342   fIntFlowCorrelationsAllPro->Fill(28.5,seven2n1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.));
2343  } // end of if(dMult>6)
2344  
2345  // 8-particle:
2346  Double_t eight1n1n1n1n1n1n1n1n = 0.; // <cos(n*(phi1+phi2+phi3+phi4-phi5-phi6-phi7-phi8))>
2347  if(dMult>7)
2348  {
2349   eight1n1n1n1n1n1n1n1n = (pow(pow(dReQ1n,2.)+pow(dImQ1n,2.),4.)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)
2350                         * (12.*seven2n1n1n1n1n1n1n+16.*six1n1n1n1n1n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)
2351                         * (8.*six3n1n1n1n1n1n+48.*six1n1n1n1n1n1n+6.*six2n2n1n1n1n1n+96.*five2n1n1n1n1n+72.*four1n1n1n1n+36.*six2n1n1n2n1n1n)
2352                         - dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(2.*five4n1n1n1n1n+32.*five2n1n1n1n1n+36.*four1n1n1n1n
2353                         + 32.*four3n1n1n1n+48.*five2n1n1n1n1n+48.*five3n1n2n1n1n+144.*five2n1n1n1n1n+288.*four1n1n1n1n+36.*five2n2n2n1n1n
2354                         + 144.*three2n1n1n+96.*two1n1n+144.*four2n1n2n1n)-dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)
2355                         * (8.*four3n1n1n1n+48.*four1n1n1n1n+12.*four4n2n1n1n+96.*four2n1n2n1n+96.*three2n1n1n+72.*three2n1n1n+144.*two1n1n
2356                         + 16.*four3n1n3n1n+48.*four3n1n1n1n+144.*four1n1n1n1n+72.*four1n1n1n1n+96.*three3n2n1n+24.*four3n1n2n2n+144.*four2n1n2n1n
2357                         + 288.*two1n1n+288.*three2n1n1n+9.*four2n2n2n2n+72.*two2n2n+24.)-dMult*(dMult-1.)*(dMult-2.)*(12.*three2n1n1n+16.*two1n1n
2358                         + 24.*three3n2n1n+48.*three2n1n1n+96.*two1n1n+8.*three4n3n1n+32.*three3n2n1n+96.*three2n1n1n+144.*two1n1n+6.*three4n2n2n
2359                         + 96.*two2n2n+36.*two2n2n+72.+48.*three3n2n1n+16.*two3n3n+72.*three2n1n1n+144.*two1n1n)-dMult*(dMult-1.)*(8.*two1n1n
2360                         + 12.*two2n2n+16.+8.*two3n3n+48.*two1n1n+1.*two4n4n+16.*two2n2n+18.)-dMult)
2361                         / (dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.)); // to be improved (direct formula needed)
2362   
2363   // average 8-particle correlations for single event: 
2364   fIntFlowCorrelationsAllEBE->SetBinContent(31,eight1n1n1n1n1n1n1n1n);
2365        
2366   // average 8-particle correlations for all events:                       
2367   fIntFlowCorrelationsAllPro->Fill(30.5,eight1n1n1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.));
2368  
2369   // store separetately <8> (to be improved: do I really need this?)
2370   fIntFlowCorrelationsEBE->SetBinContent(4,eight1n1n1n1n1n1n1n1n); // <8>
2371   fIntFlowEventWeightsForCorrelationsEBE->SetBinContent(4,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.)); // eW_<8>
2372   fIntFlowCorrelationsPro->Fill(3.5,eight1n1n1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.));
2373   
2374   // distribution of <cos(n*(phi1+phi2+phi3+phi4-phi5-phi6-phi7-phi8))>
2375   //f8pDistribution->Fill(eight1n1n1n1n1n1n1n1n,dMult*(dMult-1.)*(dMult-2.)*(dMult-3.)*(dMult-4.)*(dMult-5.)*(dMult-6.)*(dMult-7.));
2376  } // end of if(dMult>7) 
2377  
2378 } // end of AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
2379
2380
2381 //================================================================================================================================
2382
2383
2384 void AliFlowAnalysisWithQCumulants::CalculateIntFlowProductOfCorrelations()
2385 {
2386  // Calculate averages of products of correlations for integrated flow  // to be improved (this method can be implemented better)
2387  
2388  // a) Binning of fIntFlowProductOfCorrelationsPro is organized as follows:
2389  //     1st bin: <<2><4>> 
2390  //     2nd bin: <<2><6>>
2391  //     3rd bin: <<2><8>>
2392  //     4th bin: <<4><6>>
2393  //     5th bin: <<4><8>>
2394  //     6th bin: <<6><8>>
2395
2396  /*
2397  Double_t dMult = (*fSMpk)(0,0); // multiplicity 
2398
2399  Double_t twoEBE = fIntFlowCorrelationsEBE->GetBinContent(1); // <2>
2400  Double_t fourEBE = fIntFlowCorrelationsEBE->GetBinContent(2); // <4>
2401  Double_t sixEBE = fIntFlowCorrelationsEBE->GetBinContent(3); // <6>
2402  Double_t eightEBE = fIntFlowCorrelationsEBE->GetBinContent(4); // <8>
2403  
2404  Double_t eW2 = 0.; // event weight for <2>
2405  Double_t eW4 = 0.; // event weight for <4>
2406  Double_t eW6 = 0.; // event weight for <6>
2407  Double_t eW8 = 0.; // event weight for <8>
2408  
2409  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
2410  {
2411   eW2 = dMult*(dMult-1);
2412   eW4 = dMult*(dMult-1)*(dMult-2)*(dMult-3);
2413   eW6 = dMult*(dMult-1)*(dMult-2)*(dMult-3)*(dMult-4)*(dMult-5);
2414   eW8 = dMult*(dMult-1)*(dMult-2)*(dMult-3)*(dMult-4)*(dMult-5)*(dMult-6)*(dMult-7);
2415  } else 
2416    {
2417     eW2 = (*fSMpk)(1,1)-(*fSMpk)(0,2); // dM11 = sum_{i,j=1,i!=j}^M w_i w_j;
2418     eW4 = (*fSMpk)(3,1)-6.*(*fSMpk)(0,2)*(*fSMpk)(1,1)  
2419         + 8.*(*fSMpk)(0,3)*(*fSMpk)(0,1)
2420         + 3.*(*fSMpk)(1,2)-6.*(*fSMpk)(0,4); // dM1111 = sum_{i,j,k,l=1,i!=j!=k!=l}^M w_i w_j w_k w_l
2421    }
2422   
2423  fIntFlowProductOfCorrelationsPro->Fill(0.5,twoEBE*fourEBE,eW2*eW4); // <<2><4>> 
2424  fIntFlowProductOfCorrelationsPro->Fill(1.5,twoEBE*sixEBE,eW2*eW6); // <<2><6>>
2425  fIntFlowProductOfCorrelationsPro->Fill(2.5,twoEBE*eightEBE,eW2*eW8); // <<2><8>>
2426  fIntFlowProductOfCorrelationsPro->Fill(3.5,fourEBE*sixEBE,eW4*eW6); // <<4><6>>
2427  fIntFlowProductOfCorrelationsPro->Fill(4.5,fourEBE*eightEBE,eW4*eW8); // <<4><8>>
2428  fIntFlowProductOfCorrelationsPro->Fill(5.5,sixEBE*eightEBE,eW6*eW8); // <<6><8>>
2429  */
2430  
2431  
2432  Int_t counter = 0;
2433  
2434  for(Int_t ci1=1;ci1<4;ci1++)
2435  {
2436   for(Int_t ci2=ci1+1;ci2<=4;ci2++)
2437   {
2438    fIntFlowProductOfCorrelationsPro->Fill(0.5+counter++,
2439                                           fIntFlowCorrelationsEBE->GetBinContent(ci1)*fIntFlowCorrelationsEBE->GetBinContent(ci2),
2440                                           fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci1)*fIntFlowEventWeightsForCorrelationsEBE->GetBinContent(ci2));
2441   }
2442  }
2443  
2444 } // end of AliFlowAnalysisWithQCumulants::CalculateIntFlowProductOfCorrelations()
2445
2446
2447 //================================================================================================================================
2448
2449
2450 void AliFlowAnalysisWithQCumulants::CalculateCovariancesIntFlow()
2451 {
2452  // a) Calculate unbiased estimators Cov(<2>,<4>), Cov(<2>,<6>), Cov(<2>,<8>), Cov(<4>,<6>), Cov(<4>,<8>) and Cov(<6>,<8>)
2453  //    for covariances V_(<2>,<4>), V_(<2>,<6>), V_(<2>,<8>), V_(<4>,<6>), V_(<4>,<8>) and V_(<6>,<8>).
2454  // b) Store in histogram fIntFlowCovariances for instance the following: 
2455  //
2456  //             Cov(<2>,<4>) * (sum_{i=1}^{N} w_{<2>}_i w_{<4>}_i )/[(sum_{i=1}^{N} w_{<2>}_i) * (sum_{j=1}^{N} w_{<4>}_j)]
2457  // 
2458  //    where N is the number of events, w_{<2>} is event weight for <2> and w_{<4>} is event weight for <4>.
2459  // c) Binning of fIntFlowCovariances is organized as follows:
2460  // 
2461  //     1st bin: Cov(<2>,<4>) * (sum_{i=1}^{N} w_{<2>}_i w_{<4>}_i )/[(sum_{i=1}^{N} w_{<2>}_i) * (sum_{j=1}^{N} w_{<4>}_j)] 
2462  //     2nd bin: Cov(<2>,<6>) * (sum_{i=1}^{N} w_{<2>}_i w_{<6>}_i )/[(sum_{i=1}^{N} w_{<2>}_i) * (sum_{j=1}^{N} w_{<6>}_j)]
2463  //     3rd bin: Cov(<2>,<8>) * (sum_{i=1}^{N} w_{<2>}_i w_{<8>}_i )/[(sum_{i=1}^{N} w_{<2>}_i) * (sum_{j=1}^{N} w_{<8>}_j)]
2464  //     4th bin: Cov(<4>,<6>) * (sum_{i=1}^{N} w_{<4>}_i w_{<6>}_i )/[(sum_{i=1}^{N} w_{<4>}_i) * (sum_{j=1}^{N} w_{<6>}_j)]
2465  //     5th bin: Cov(<4>,<8>) * (sum_{i=1}^{N} w_{<4>}_i w_{<8>}_i )/[(sum_{i=1}^{N} w_{<4>}_i) * (sum_{j=1}^{N} w_{<8>}_j)]
2466  //     6th bin: Cov(<6>,<8>) * (sum_{i=1}^{N} w_{<6>}_i w_{<8>}_i )/[(sum_{i=1}^{N} w_{<6>}_i) * (sum_{j=1}^{N} w_{<8>}_j)]
2467     
2468  for(Int_t power=0;power<2;power++)
2469  { 
2470   if(!(fIntFlowCorrelationsPro && fIntFlowProductOfCorrelationsPro 
2471        && fIntFlowSumOfEventWeights[power] && fIntFlowSumOfProductOfEventWeights
2472        && fIntFlowCovariances)) 
2473   {
2474    cout<<"WARNING: fIntFlowCorrelationsPro && fIntFlowProductOfCorrelationsPro "<<endl;
2475    cout<<"         && fIntFlowSumOfEventWeights[power] && fIntFlowSumOfProductOfEventWeights"<<endl;
2476    cout<<"         && fIntFlowCovariances is NULL in AFAWQC::FCIF() !!!!"<<endl;
2477    cout<<"power = "<<power<<endl;
2478    exit(0);
2479   }
2480  }
2481    
2482  // average 2-, 4-, 6- and 8-particle correlations for all events:
2483  Double_t correlation[4] = {0.};
2484  for(Int_t ci=0;ci<4;ci++)
2485  {
2486   correlation[ci] = fIntFlowCorrelationsPro->GetBinContent(ci+1);
2487  } 
2488  // average products of 2-, 4-, 6- and 8-particle correlations: 
2489  Double_t productOfCorrelations[4][4] = {{0.}};
2490  Int_t productOfCorrelationsLabel = 1;
2491  // denominators in the expressions for the unbiased estimator for covariance:
2492  Double_t denominator[4][4] = {{0.}};
2493  Int_t sumOfProductOfEventWeightsLabel1 = 1;
2494  // weight dependent prefactor which multiply unbiased estimators for covariances:
2495  Double_t wPrefactor[4][4] = {{0.}}; 
2496  Int_t sumOfProductOfEventWeightsLabel2 = 1;
2497  for(Int_t c1=0;c1<4;c1++)
2498  {
2499   for(Int_t c2=c1+1;c2<4;c2++)
2500   {
2501    productOfCorrelations[c1][c2] = fIntFlowProductOfCorrelationsPro->GetBinContent(productOfCorrelationsLabel);
2502    if(fIntFlowSumOfEventWeights[0]->GetBinContent(c1+1) && fIntFlowSumOfEventWeights[0]->GetBinContent(c2+1))
2503    {
2504     denominator[c1][c2] = 1.-(fIntFlowSumOfProductOfEventWeights->GetBinContent(sumOfProductOfEventWeightsLabel1))/
2505                              (fIntFlowSumOfEventWeights[0]->GetBinContent(c1+1) 
2506                               * fIntFlowSumOfEventWeights[0]->GetBinContent(c2+1));
2507                               
2508     wPrefactor[c1][c2] =  fIntFlowSumOfProductOfEventWeights->GetBinContent(sumOfProductOfEventWeightsLabel2)/ 
2509                           (fIntFlowSumOfEventWeights[0]->GetBinContent(c1+1)
2510                             * fIntFlowSumOfEventWeights[0]->GetBinContent(c2+1));
2511                           
2512                               
2513    }
2514    productOfCorrelationsLabel++;
2515    sumOfProductOfEventWeightsLabel1++;
2516    sumOfProductOfEventWeightsLabel2++;  
2517   }
2518  }
2519  
2520  // covariance label:
2521  Int_t covarianceLabel = 1;
2522  for(Int_t c1=0;c1<4;c1++)
2523  {
2524   for(Int_t c2=c1+1;c2<4;c2++)
2525   {
2526    if(denominator[c1][c2])
2527    {
2528     // covariances:
2529     Double_t cov = (productOfCorrelations[c1][c2]-correlation[c1]*correlation[c2])/denominator[c1][c2]; 
2530     // covarianced multiplied with weight dependent prefactor:
2531     Double_t wCov = cov * wPrefactor[c1][c2];
2532     fIntFlowCovariances->SetBinContent(covarianceLabel,wCov);
2533    }
2534    covarianceLabel++;
2535   }
2536  }
2537  
2538 } // end of AliFlowAnalysisWithQCumulants::CalculateCovariancesIntFlow()
2539
2540
2541 //================================================================================================================================
2542
2543
2544 void AliFlowAnalysisWithQCumulants::FinalizeCorrelationsIntFlow() 
2545 {
2546  // From profile fIntFlowCorrelationsPro access measured correlations and spread, 
2547  // correctly calculate the statistical errors and store the final results and 
2548  // statistical errors for correlations in histogram fIntFlowCorrelationsHist.
2549  //
2550  // Remark: Statistical error of correlation is calculated as:
2551  //
2552  //          statistical error = termA * spread * termB:
2553  //          termA = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
2554  //          termB = 1/sqrt(1-termA^2)   
2555  
2556  for(Int_t power=0;power<2;power++)
2557  { 
2558   if(!(fIntFlowCorrelationsHist && fIntFlowCorrelationsPro && fIntFlowSumOfEventWeights[power])) 
2559   {
2560    cout<<"WARNING: fIntFlowCorrelationsHist && fIntFlowCorrelationsPro && fIntFlowSumOfEventWeights[power] is NULL in AFAWQC::FCIF() !!!!"<<endl;
2561    cout<<"power = "<<power<<endl;
2562    exit(0);
2563   }
2564  }
2565   
2566  for(Int_t ci=1;ci<=4;ci++) // correlation index
2567  {
2568   Double_t correlation = fIntFlowCorrelationsPro->GetBinContent(ci);
2569   Double_t spread = fIntFlowCorrelationsPro->GetBinError(ci);
2570   Double_t sumOfLinearEventWeights = fIntFlowSumOfEventWeights[0]->GetBinContent(ci);
2571   Double_t sumOfQuadraticEventWeights = fIntFlowSumOfEventWeights[1]->GetBinContent(ci);
2572   Double_t termA = 0.;
2573   Double_t termB = 0.;
2574   if(sumOfLinearEventWeights)
2575   {
2576    termA = pow(sumOfQuadraticEventWeights,0.5)/sumOfLinearEventWeights;
2577   } else
2578     {
2579      cout<<"WARNING: sumOfLinearEventWeights == 0 in AFAWQC::FCIF() !!!!"<<endl;
2580      cout<<"         (for "<<2*ci<<"-particle correlation)"<<endl;
2581     }
2582   if(1.-pow(termA,2.) > 0.)
2583   {
2584    termB = 1./pow(1-pow(termA,2.),0.5);
2585   } else
2586     {
2587      cout<<"WARNING: 1.-pow(termA,2.) <= 0 in AFAWQC::FCIF() !!!!"<<endl;   
2588      cout<<"         (for "<<2*ci<<"-particle correlation)"<<endl;
2589     }     
2590   Double_t statisticalError = termA * spread * termB;
2591   fIntFlowCorrelationsHist->SetBinContent(ci,correlation);
2592   fIntFlowCorrelationsHist->SetBinError(ci,statisticalError);
2593  } // end of for(Int_t ci=1;ci<=4;ci++) // correlation index                                                                
2594                                                                                                                               
2595 } // end of AliFlowAnalysisWithQCumulants::FinalizeCorrelationsIntFlow()
2596
2597
2598 //================================================================================================================================
2599
2600
2601 void AliFlowAnalysisWithQCumulants::FillAverageMultiplicities(Int_t nRP)
2602 {
2603  // Fill profile fAverageMultiplicity to hold average multiplicities and number of events for events with nRP>=0, nRP>=1, ... , and nRP>=8
2604  
2605  // Binning of fAverageMultiplicity is organized as follows:
2606  //  1st bin: all events (including the empty ones)
2607  //  2nd bin: event with # of RPs greater or equal to 1
2608  //  3rd bin: event with # of RPs greater or equal to 2
2609  //  4th bin: event with # of RPs greater or equal to 3
2610  //  5th bin: event with # of RPs greater or equal to 4
2611  //  6th bin: event with # of RPs greater or equal to 5
2612  //  7th bin: event with # of RPs greater or equal to 6
2613  //  8th bin: event with # of RPs greater or equal to 7
2614  //  9th bin: event with # of RPs greater or equal to 8
2615  
2616  if(!fAvMultiplicity)
2617  {
2618   cout<<"WARNING: fAvMultiplicity is NULL in AFAWQC::FAM() !!!!"<<endl;
2619   exit(0);
2620  }
2621  
2622  if(nRP<0)
2623  {
2624   cout<<"WARNING: nRP<0 in in AFAWQC::FAM() !!!!"<<endl;
2625   exit(0);
2626  }
2627  
2628  for(Int_t i=0;i<9;i++)
2629  {
2630   if(nRP>=i) fAvMultiplicity->Fill(i+0.5,nRP,1);
2631  }
2632  
2633 } // end of AliFlowAnalysisWithQCumulants::FillAverageMultiplicities(nRP)
2634
2635
2636 //================================================================================================================================
2637
2638
2639 void AliFlowAnalysisWithQCumulants::CalculateCumulantsIntFlow()
2640 {
2641  // a) Calculate Q-cumulants from the measured multiparticle correlations.
2642  // b) Propagate the statistical errors of measured multiparticle correlations to statistical errors of Q-cumulants.  
2643  // c) REMARK: Q-cumulants calculated in this method are biased by non-uniform acceptance of detector !!!! 
2644  //            Method ApplyCorrectionForNonUniformAcceptance* (to be improved: finalize the name here)
2645  //            is called afterwards to correct for this bias.   
2646  // d) Store the results and statistical error of Q-cumulants in histogram fCumulants.
2647  //    Binning of fCumulants is organized as follows:
2648  //
2649  //     1st bin: QC{2}
2650  //     2nd bin: QC{4}
2651  //     3rd bin: QC{6}
2652  //     4th bin: QC{8}
2653  
2654  if(!(fIntFlowCorrelationsHist && fIntFlowCovariances && fIntFlowQcumulants))
2655  {
2656   cout<<"WARNING: fIntFlowCorrelationsHist && fIntFlowCovariances && fIntFlowQcumulants is NULL in AFAWQC::CCIF() !!!!"<<endl;
2657   exit(0);
2658  }
2659  
2660  // correlations:
2661  Double_t two = fIntFlowCorrelationsHist->GetBinContent(1); // <<2>> 
2662  Double_t four = fIntFlowCorrelationsHist->GetBinContent(2); // <<4>>  
2663  Double_t six = fIntFlowCorrelationsHist->GetBinContent(3); // <<6>> 
2664  Double_t eight = fIntFlowCorrelationsHist->GetBinContent(4); // <<8>>  
2665  
2666  // statistical errors of average 2-, 4-, 6- and 8-particle azimuthal correlations:
2667  Double_t twoError = fIntFlowCorrelationsHist->GetBinError(1); // statistical error of <2>  
2668  Double_t fourError = fIntFlowCorrelationsHist->GetBinError(2); // statistical error of <4>   
2669  Double_t sixError = fIntFlowCorrelationsHist->GetBinError(3); // statistical error of <6> 
2670  Double_t eightError = fIntFlowCorrelationsHist->GetBinError(4); // statistical error of <8> 
2671  
2672  // covariances (multiplied by prefactor depending on weights - see comments in CalculateCovariancesIntFlow()):
2673  Double_t wCov24 = fIntFlowCovariances->GetBinContent(1); // Cov(<2>,<4>) * prefactor(w_<2>,w_<4>)
2674  Double_t wCov26 = fIntFlowCovariances->GetBinContent(2); // Cov(<2>,<6>) * prefactor(w_<2>,w_<6>)
2675  Double_t wCov28 = fIntFlowCovariances->GetBinContent(3); // Cov(<2>,<8>) * prefactor(w_<2>,w_<8>)
2676  Double_t wCov46 = fIntFlowCovariances->GetBinContent(4); // Cov(<4>,<6>) * prefactor(w_<4>,w_<6>)
2677  Double_t wCov48 = fIntFlowCovariances->GetBinContent(5); // Cov(<4>,<8>) * prefactor(w_<4>,w_<8>)
2678  Double_t wCov68 = fIntFlowCovariances->GetBinContent(6); // Cov(<6>,<8>) * prefactor(w_<6>,w_<8>)
2679  
2680  // Q-cumulants: 
2681  Double_t qc2 = 0.; // QC{2}
2682  Double_t qc4 = 0.; // QC{4}
2683  Double_t qc6 = 0.; // QC{6}
2684  Double_t qc8 = 0.; // QC{8}
2685  if(two) qc2 = two; 
2686  if(four) qc4 = four-2.*pow(two,2.); 
2687  if(six) qc6 = six-9.*two*four+12.*pow(two,3.); 
2688  if(eight) qc8 = eight-16.*two*six-18.*pow(four,2.)+144.*pow(two,2.)*four-144.*pow(two,4.); 
2689  
2690  // statistical errors of Q-cumulants:       
2691  Double_t qc2Error = 0.;
2692  Double_t qc4Error = 0.;
2693  Double_t qc6Error = 0.;
2694  Double_t qc8Error = 0.;
2695  
2696  // squared statistical errors of Q-cumulants:       
2697  //Double_t qc2ErrorSquared = 0.;
2698  Double_t qc4ErrorSquared = 0.;
2699  Double_t qc6ErrorSquared = 0.;
2700  Double_t qc8ErrorSquared = 0.;
2701         
2702  // statistical error of QC{2}:              
2703  qc2Error = twoError;                     
2704                              
2705  // statistical error of QC{4}:              
2706  qc4ErrorSquared = 16.*pow(two,2.)*pow(twoError,2)+pow(fourError,2.)
2707                  - 8.*two*wCov24;                     
2708  if(qc4ErrorSquared>0.)
2709  {
2710   qc4Error = pow(qc4ErrorSquared,0.5);
2711  } else 
2712    {
2713     cout<<"WARNING: Statistical error of QC{4} is imaginary !!!!"<<endl;
2714    }
2715                                            
2716  // statistical error of QC{6}:              
2717  qc6ErrorSquared = 81.*pow(4.*pow(two,2.)-four,2.)*pow(twoError,2.)
2718                  + 81.*pow(two,2.)*pow(fourError,2.)
2719                  + pow(sixError,2.)
2720                  - 162.*two*(4.*pow(two,2.)-four)*wCov24
2721                  + 18.*(4.*pow(two,2.)-four)*wCov26
2722                  - 18.*two*wCov46; 
2723                     
2724  if(qc6ErrorSquared>0.)
2725  {
2726   qc6Error = pow(qc6ErrorSquared,0.5);
2727  } else 
2728    {
2729     cout<<"WARNING: Statistical error of QC{6} is imaginary !!!!"<<endl;
2730    }
2731                             
2732  // statistical error of QC{8}:              
2733  qc8ErrorSquared = 256.*pow(36.*pow(two,3.)-18.*four*two+six,2.)*pow(twoError,2.)
2734                  + 1296.*pow(4.*pow(two,2.)-four,2.)*pow(fourError,2.)
2735                  + 256.*pow(two,2.)*pow(sixError,2.)
2736                  + pow(eightError,2.)
2737                  - 1152.*(36.*pow(two,3.)-18.*four*two+six)*(4.*pow(two,2.)-four)*wCov24
2738                  + 512.*two*(36.*pow(two,3.)-18.*four*two+six)*wCov26
2739                  - 32.*(36.*pow(two,3.)-18.*four*two+six)*wCov28
2740                  - 1152.*two*(4.*pow(two,2.)-four)*wCov46
2741                  + 72.*(4.*pow(two,2.)-four)*wCov48
2742                  - 32.*two*wCov68;      
2743  if(qc8ErrorSquared>0.)
2744  {
2745   qc8Error = pow(qc8ErrorSquared,0.5);
2746  } else 
2747    {
2748     cout<<"WARNING: Statistical error of QC{8} is imaginary !!!!"<<endl;
2749    }
2750
2751  // store the results and statistical errors for Q-cumulants:
2752  fIntFlowQcumulants->SetBinContent(1,qc2);
2753  fIntFlowQcumulants->SetBinError(1,qc2Error);
2754  fIntFlowQcumulants->SetBinContent(2,qc4);
2755  fIntFlowQcumulants->SetBinError(2,qc4Error);
2756  fIntFlowQcumulants->SetBinContent(3,qc6);
2757  fIntFlowQcumulants->SetBinError(3,qc6Error);
2758  fIntFlowQcumulants->SetBinContent(4,qc8); 
2759  fIntFlowQcumulants->SetBinError(4,qc8Error); 
2760   
2761 } // end of AliFlowAnalysisWithQCumulants::CalculateCumulantsIntFlow()
2762
2763
2764 //================================================================================================================================ 
2765
2766
2767 void AliFlowAnalysisWithQCumulants::CalculateIntFlow()
2768 {
2769  // a) Calculate the final results for integrated flow estimates from Q-cumulants.
2770  // b) Propagate the statistical errors of measured multiparticle correlations to statistical errors of integrated flow estimates.  
2771  // c) Store the results and statistical errors of integrated flow estimates in histogram fIntFlow.
2772  //    Binning of fIntFlow is organized as follows:
2773  //
2774  //     1st bin: v{2,QC}
2775  //     2nd bin: v{4,QC}
2776  //     3rd bin: v{6,QC}
2777  //     4th bin: v{8,QC}
2778  
2779  if(!(fIntFlowCorrelationsHist && fIntFlowCovariances && fIntFlowQcumulants && fIntFlow))
2780  {
2781   cout<<"WARNING: fIntFlowCorrelationsHist && fIntFlowCovariances && fIntFlowQcumulants && fIntFlow is NULL in AFAWQC::CCIF() !!!!"<<endl;
2782   exit(0);
2783  }
2784    
2785  // Q-cumulants:
2786  Double_t qc2 = fIntFlowQcumulants->GetBinContent(1); // QC{2}  
2787  Double_t qc4 = fIntFlowQcumulants->GetBinContent(2); // QC{4}  
2788  Double_t qc6 = fIntFlowQcumulants->GetBinContent(3); // QC{6}  
2789  Double_t qc8 = fIntFlowQcumulants->GetBinContent(4); // QC{8}
2790   
2791  // correlations:
2792  Double_t two = fIntFlowCorrelationsHist->GetBinContent(1); // <<2>> 
2793  Double_t four = fIntFlowCorrelationsHist->GetBinContent(2); // <<4>>  
2794  Double_t six = fIntFlowCorrelationsHist->GetBinContent(3); // <<6>> 
2795  Double_t eight = fIntFlowCorrelationsHist->GetBinContent(4); // <<8>>  
2796  
2797  // statistical errors of average 2-, 4-, 6- and 8-particle azimuthal correlations:
2798  Double_t twoError = fIntFlowCorrelationsHist->GetBinError(1); // statistical error of <2>  
2799  Double_t fourError = fIntFlowCorrelationsHist->GetBinError(2); // statistical error of <4>   
2800  Double_t sixError = fIntFlowCorrelationsHist->GetBinError(3); // statistical error of <6> 
2801  Double_t eightError = fIntFlowCorrelationsHist->GetBinError(4); // statistical error of <8> 
2802  
2803  // covariances (multiplied by prefactor depending on weights - see comments in CalculateCovariancesIntFlow()):
2804  Double_t wCov24 = fIntFlowCovariances->GetBinContent(1); // Cov(<2>,<4>) * prefactor(w_<2>,w_<4>)
2805  Double_t wCov26 = fIntFlowCovariances->GetBinContent(2); // Cov(<2>,<6>) * prefactor(w_<2>,w_<6>)
2806  Double_t wCov28 = fIntFlowCovariances->GetBinContent(3); // Cov(<2>,<8>) * prefactor(w_<2>,w_<8>)
2807  Double_t wCov46 = fIntFlowCovariances->GetBinContent(4); // Cov(<4>,<6>) * prefactor(w_<4>,w_<6>)
2808  Double_t wCov48 = fIntFlowCovariances->GetBinContent(5); // Cov(<4>,<8>) * prefactor(w_<4>,w_<8>)
2809  Double_t wCov68 = fIntFlowCovariances->GetBinContent(6); // Cov(<6>,<8>) * prefactor(w_<6>,w_<8>)
2810   
2811  // integrated flow estimates:
2812  Double_t v2 = 0.; // v{2,QC}  
2813  Double_t v4 = 0.; // v{4,QC}  
2814  Double_t v6 = 0.; // v{6,QC}  
2815  Double_t v8 = 0.; // v{8,QC}
2816  
2817  // calculate integrated flow estimates from Q-cumulants: 
2818  if(qc2>=0.) v2 = pow(qc2,1./2.); 
2819  if(qc4<=0.) v4 = pow(-1.*qc4,1./4.); 
2820  if(qc6>=0.) v6 = pow((1./4.)*qc6,1./6.); 
2821  if(qc8<=0.) v8 = pow((-1./33.)*qc8,1./8.); 
2822    
2823  // statistical errors of integrated flow estimates:
2824  Double_t v2Error = 0.; // statistical error of v{2,QC}  
2825  Double_t v4Error = 0.; // statistical error of v{4,QC}  
2826  Double_t v6Error = 0.; // statistical error of v{6,QC}  
2827  Double_t v8Error = 0.; // statistical error of v{8,QC}
2828    
2829  // squares of statistical errors of integrated flow estimates:
2830  Double_t v2ErrorSquared = 0.; // squared statistical error of v{2,QC} 
2831  Double_t v4ErrorSquared = 0.; // squared statistical error of v{4,QC}   
2832  Double_t v6ErrorSquared = 0.; // squared statistical error of v{6,QC}   
2833  Double_t v8ErrorSquared = 0.; // squared statistical error of v{8,QC} 
2834  
2835  // calculate squared statistical errors of integrated flow estimates:
2836  if(two != 0.) 
2837  { 
2838   v2ErrorSquared = (1./(4.*two))*pow(twoError,2.);
2839  } 
2840  if(2.*pow(two,2.)-four > 0.)
2841  {
2842   v4ErrorSquared = (1./pow(2.*pow(two,2.)-four,3./2.))*
2843                    (pow(two,2.)*pow(twoError,2.)+(1./16.)*pow(fourError,2.)-(1./2.)*two*wCov24);
2844  }
2845  if(six-9.*four*two+12.*pow(two,3.) > 0.) 
2846  {
2847   v6ErrorSquared = ((1./2.)*(1./pow(2.,2./3.))*(1./pow(six-9.*four*two+12.*pow(two,3.),5./3.)))*
2848                    ((9./2.)*pow(4.*pow(two,2.)-four,2.)*pow(twoError,2.) 
2849                     + (9./2.)*pow(two,2.)*pow(fourError,2.)+(1./18.)*pow(sixError,2.)
2850                     - 9.*two*(4.*pow(two,2.)-four)*wCov24+(4.*pow(two,2.)-four)*wCov26-two*wCov46); 
2851  }
2852  if(-1.*eight+16.*six*two+18.*pow(four,2.)-144.*four*pow(two,2.)+144.*pow(two,4.) > 0.) 
2853  {
2854   v8ErrorSquared = (4./pow(33,1./4.))*(1./pow(-1.*eight+16.*six*two+18.*pow(four,2.)-144.*four*pow(two,2.)+144.*pow(two,4.),7./4.))*
2855                    (pow(36.*pow(two,3.)-18.*four*two+six,2.)*pow(twoError,2.)
2856                     + (81./16.)*pow(4.*pow(two,2.)-four,2.)*pow(fourError,2.)
2857                     + pow(two,2.)*pow(sixError,2.)
2858                     + (1./256.)*pow(eightError,2.)
2859                     - (9./2.)*(36.*pow(two,3.)-18.*four*two+six)*(4.*pow(two,2.)-four)*wCov24
2860                     + 2.*two*(36.*pow(two,3.)-18.*four*two+six)*wCov26
2861                     - (1./8.)*(36.*pow(two,3.)-18.*four*two+six)*wCov28                    
2862                     - (9./2.)*two*(4.*pow(two,2.)-four)*wCov46                   
2863                     + (9./32.)*(4.*pow(two,2.)-four)*wCov48                    
2864                     - (1./8.)*two*wCov68);
2865  } 
2866
2867  // calculate statistical errors of integrated flow estimates: 
2868  if(v2ErrorSquared > 0.)
2869  {
2870   v2Error = pow(v2ErrorSquared,0.5);
2871  } else
2872    {
2873     cout<<"WARNING: Statistical error of v{2,QC} is imaginary !!!!"<<endl;
2874    }    
2875  if(v4ErrorSquared > 0.)
2876  {
2877   v4Error = pow(v4ErrorSquared,0.5);
2878  } else
2879    {
2880     cout<<"WARNING: Statistical error of v{4,QC} is imaginary !!!!"<<endl;
2881    }     
2882  if(v6ErrorSquared > 0.)
2883  {
2884   v6Error = pow(v6ErrorSquared,0.5);
2885  } else
2886    {
2887     cout<<"WARNING: Statistical error of v{6,QC} is imaginary !!!!"<<endl;
2888    }     
2889  if(v8ErrorSquared > 0.)
2890  {
2891   v8Error = pow(v8ErrorSquared,0.5);
2892  } else
2893    {
2894     cout<<"WARNING: Statistical error of v{8,QC} is imaginary !!!!"<<endl;
2895    }    
2896                      
2897  // store the results and statistical errors of integrated flow estimates:
2898  fIntFlow->SetBinContent(1,v2);
2899  fIntFlow->SetBinError(1,v2Error);
2900  fIntFlow->SetBinContent(2,v4);
2901  fIntFlow->SetBinError(2,v4Error);
2902  fIntFlow->SetBinContent(3,v6);
2903  fIntFlow->SetBinError(3,v6Error);
2904  fIntFlow->SetBinContent(4,v8);
2905  fIntFlow->SetBinError(4,v8Error);
2906        
2907 } // end of AliFlowAnalysisWithQCumulants::CalculateIntFlow()
2908
2909
2910 //================================================================================================================================ 
2911
2912
2913 void AliFlowAnalysisWithQCumulants::FillCommonHistResultsIntFlow()
2914 {
2915  // Fill in AliFlowCommonHistResults histograms relevant for 'NONAME' integrated flow (to be improved (name))
2916  
2917  if(!fIntFlow)
2918  {
2919   cout<<"WARNING: fIntFlow is NULL in AFAWQC::FCHRIF() !!!!"<<endl;
2920   exit(0); 
2921  }  
2922     
2923  if(!(fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && fCommonHistsResults8th))
2924  {
2925   cout<<"WARNING: fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && fCommonHistsResults8th"<<endl; 
2926   cout<<"         is NULL in AFAWQC::FCHRIF() !!!!"<<endl;
2927   exit(0);
2928  }
2929  
2930  Double_t v2 = fIntFlow->GetBinContent(1);
2931  Double_t v4 = fIntFlow->GetBinContent(2);
2932  Double_t v6 = fIntFlow->GetBinContent(3);
2933  Double_t v8 = fIntFlow->GetBinContent(4);
2934   
2935  Double_t v2Error = fIntFlow->GetBinError(1);
2936  Double_t v4Error = fIntFlow->GetBinError(2);
2937  Double_t v6Error = fIntFlow->GetBinError(3);
2938  Double_t v8Error = fIntFlow->GetBinError(4);
2939  
2940  fCommonHistsResults2nd->FillIntegratedFlow(v2,v2Error); // to be improved (hardwired 2nd in the name)  
2941  fCommonHistsResults4th->FillIntegratedFlow(v4,v4Error); // to be improved (hardwired 4th in the name)
2942  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)) // to be improved (calculate also 6th and 8th order)
2943  {
2944   fCommonHistsResults6th->FillIntegratedFlow(v6,v6Error); // to be improved (hardwired 6th in the name)
2945   fCommonHistsResults8th->FillIntegratedFlow(v8,v8Error); // to be improved (hardwired 8th in the name) 
2946  }
2947  
2948 } // end of AliFlowAnalysisWithQCumulants::FillCommonHistResultsIntFlow()
2949
2950
2951 //================================================================================================================================ 
2952
2953
2954 /*
2955 void AliFlowAnalysisWithQCumulants::ApplyCorrectionForNonUniformAcceptanceToCumulantsForIntFlow(Bool_t useParticleWeights, TString eventWeights)
2956 {
2957  // apply correction for non-uniform acceptance to cumulants for integrated flow 
2958  // (Remark: non-corrected cumulants are accessed from fCumulants[pW][0], corrected cumulants are stored in fCumulants[pW][1])
2959  
2960  // shortcuts for the flags:
2961  Int_t pW = (Int_t)(useParticleWeights); // 0=pWeights not used, 1=pWeights used
2962  Int_t eW = -1;
2963  
2964  if(eventWeights == "exact")
2965  {
2966   eW = 0;
2967  }
2968  
2969  if(!(fCumulants[pW][eW][0] && fCumulants[pW][eW][1] && fCorrections[pW][eW]))
2970  {
2971   cout<<"WARNING: fCumulants[pW][eW][0] && fCumulants[pW][eW][1] && fCorrections[pW][eW] is NULL in AFAWQC::ACFNUATCFIF() !!!!"<<endl;
2972   cout<<"pW = "<<pW<<endl;
2973   cout<<"eW = "<<eW<<endl;
2974   exit(0);
2975  } 
2976   
2977  // non-corrected cumulants:
2978  Double_t qc2 = fCumulants[pW][eW][0]->GetBinContent(1); 
2979  Double_t qc4 = fCumulants[pW][eW][0]->GetBinContent(2); 
2980  Double_t qc6 = fCumulants[pW][eW][0]->GetBinContent(3); 
2981  Double_t qc8 = fCumulants[pW][eW][0]->GetBinContent(4); 
2982  // statistical error of non-corrected cumulants:  
2983  Double_t qc2Error = fCumulants[pW][eW][0]->GetBinError(1); 
2984  Double_t qc4Error = fCumulants[pW][eW][0]->GetBinError(2); 
2985  Double_t qc6Error = fCumulants[pW][eW][0]->GetBinError(3); 
2986  Double_t qc8Error = fCumulants[pW][eW][0]->GetBinError(4); 
2987  // corrections for non-uniform acceptance:
2988  Double_t qc2Correction = fCorrections[pW][eW]->GetBinContent(1); 
2989  Double_t qc4Correction = fCorrections[pW][eW]->GetBinContent(2); 
2990  Double_t qc6Correction = fCorrections[pW][eW]->GetBinContent(3); 
2991  Double_t qc8Correction = fCorrections[pW][eW]->GetBinContent(4); 
2992  // corrected cumulants:
2993  Double_t qc2Corrected = qc2 + qc2Correction;
2994  Double_t qc4Corrected = qc4 + qc4Correction;
2995  Double_t qc6Corrected = qc6 + qc6Correction;
2996  Double_t qc8Corrected = qc8 + qc8Correction;
2997   
2998  // ... to be improved (I need here also to correct error of QCs for NUA. 
2999  // For simplicity sake I assume at the moment that this correction is negliglible, but it will be added eventually...)
3000  
3001  // store corrected results and statistical errors for cumulants:   
3002  fCumulants[pW][eW][1]->SetBinContent(1,qc2Corrected);
3003  fCumulants[pW][eW][1]->SetBinContent(2,qc4Corrected);
3004  fCumulants[pW][eW][1]->SetBinContent(3,qc6Corrected);
3005  fCumulants[pW][eW][1]->SetBinContent(4,qc8Corrected);
3006  fCumulants[pW][eW][1]->SetBinError(1,qc2Error); // to be improved (correct also qc2Error for NUA)
3007  fCumulants[pW][eW][1]->SetBinError(2,qc4Error); // to be improved (correct also qc4Error for NUA)
3008  fCumulants[pW][eW][1]->SetBinError(3,qc6Error); // to be improved (correct also qc6Error for NUA)
3009  fCumulants[pW][eW][1]->SetBinError(4,qc8Error); // to be improved (correct also qc8Error for NUA)  
3010   
3011 } // end of AliFlowAnalysisWithQCumulants::ApplyCorrectionForNonUniformAcceptanceToCumulantsForIntFlow(Bool_t useParticleWeights, TString eventWeights)
3012 */
3013
3014
3015 //================================================================================================================================
3016
3017
3018 /*  
3019 void AliFlowAnalysisWithQCumulants::PrintQuantifyingCorrectionsForNonUniformAcceptance(Bool_t useParticleWeights, TString eventWeights)
3020 {
3021  // print on the screen QC{n,biased}/QC{n,corrected}
3022  
3023  // shortcuts for the flags:
3024  Int_t pW = (Int_t)(useParticleWeights); // 0=pWeights not used, 1=pWeights used
3025  
3026  Int_t eW = -1;
3027  
3028  if(eventWeights == "exact")
3029  {
3030   eW = 0;
3031  } 
3032