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