correlators galore
[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 && fCalculateDiffFlow)
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(),58,0,58);
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",58,0,58);
1684  fIntFlowCorrelationsAllPro->Sumw2();
1685  fIntFlowCorrelationsAllPro->SetTickLength(-0.01,"Y");
1686  fIntFlowCorrelationsAllPro->SetMarkerStyle(25);
1687  fIntFlowCorrelationsAllPro->SetLabelSize(0.03);
1688  fIntFlowCorrelationsAllPro->SetLabelOffset(0.01,"Y");
1689  // 2-p correlations:
1690  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
1691  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
1692  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
1693  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
1694  // 3-p correlations:
1695  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
1696  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
1697  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
1698  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
1699  // 4-p correlations:
1700  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
1701  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
1702  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
1703  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
1704  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
1705  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
1706  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
1707  // 5-p correlations:
1708  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
1709  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
1710  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
1711  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
1712  // 6-p correlations:
1713  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
1714  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
1715  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
1716  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
1717  // 7-p correlations:  
1718  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
1719  // 8-p correlations:
1720  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
1721  //  EXTRA correlations for v3{5} study:
1722  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(33,"<<4>>_{4n,2n|3n,3n}");
1723  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(34,"<<5>>_{3n,3n|2n,2n,2n}");
1724  //  EXTRA correlations for Teaney-Yan study:
1725  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(35,"<2>_{5n|5n}");
1726  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(36,"<2>_{6n|6n}");
1727  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(37,"<3>_{5n|3n,2n}");
1728  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(38,"<3>_{5n|4n,1n}");
1729  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(39,"<3>_{6n|3n,3n}");
1730  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(40,"<3>_{6n|4n,2n}");
1731  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(41,"<3>_{6n|5n,1n}");
1732  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(42,"<4>_{6n|3n,2n,1n}");
1733  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(43,"<4>_{3n,2n|3n,2n}");
1734  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(44,"<4>_{4n,1n|3n,2n}");
1735  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(45,"<4>_{3n,3n|3n,3n}");
1736  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(46,"<4>_{4n,2n|3n,3n}");
1737  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(47,"<4>_{5n,1n|3n,3n}");
1738  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(48,"<4>_{4n,2n|4n,2n}");
1739  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(49,"<4>_{5n,1n|4n,2n}");
1740  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(50,"<4>_{5n|3n,1n,1n}");
1741  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(51,"<4>_{5n|2n,2n,1n}");
1742  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(52,"<4>_{5n,1n|5n,1n}");
1743  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(53,"<5>_{3n,3n|3n,2n,1n}");
1744  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(54,"<5>_{4n,2n|3n,2n,1n}");
1745  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(55,"<5>_{3n,2n|3n,1n,1n}");
1746  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(56,"<5>_{3n,2n|2n,2n,1n}");
1747  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(57,"<5>_{5n,1n|3n,2n,1n}");
1748  (fIntFlowCorrelationsAllPro->GetXaxis())->SetBinLabel(58,"<6>_{3n,2n,1n|3n,2n,1n}");
1749  fIntFlowProfiles->Add(fIntFlowCorrelationsAllPro);
1750  // average all correlations versus multiplicity (errors via Sumw2 - to be improved):
1751  if(fCalculateAllCorrelationsVsM)
1752  {
1753   // 2-p correlations vs M:  
1754   fIntFlowCorrelationsAllVsMPro[0] = new TProfile("two1n1n","#LT#LT2#GT#GT_{n|n}",fnBinsMult,fMinMult,fMaxMult);
1755   fIntFlowCorrelationsAllVsMPro[0]->Sumw2();
1756   fIntFlowCorrelationsAllVsMPro[0]->GetXaxis()->SetTitle("M");  
1757   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[0]);  
1758   fIntFlowCorrelationsAllVsMPro[1] = new TProfile("two2n2n","#LT#LT2#GT#GT_{2n|2n}",fnBinsMult,fMinMult,fMaxMult);
1759   fIntFlowCorrelationsAllVsMPro[1]->Sumw2();
1760   fIntFlowCorrelationsAllVsMPro[1]->GetXaxis()->SetTitle("M");  
1761   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[1]);
1762   fIntFlowCorrelationsAllVsMPro[2] = new TProfile("two3n3n","#LT#LT2#GT#GT_{3n|3n}",fnBinsMult,fMinMult,fMaxMult);
1763   fIntFlowCorrelationsAllVsMPro[2]->Sumw2();
1764   fIntFlowCorrelationsAllVsMPro[2]->GetXaxis()->SetTitle("M");  
1765   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[2]);
1766   fIntFlowCorrelationsAllVsMPro[3] = new TProfile("two4n4n","#LT#LT2#GT#GT_{4n|4n}",fnBinsMult,fMinMult,fMaxMult);
1767   fIntFlowCorrelationsAllVsMPro[3]->Sumw2();
1768   fIntFlowCorrelationsAllVsMPro[3]->GetXaxis()->SetTitle("M");  
1769   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[3]);
1770   // 3-p correlations vs M:
1771   fIntFlowCorrelationsAllVsMPro[5] = new TProfile("three2n1n1n","#LT#LT3#GT#GT_{2n|n,n}",fnBinsMult,fMinMult,fMaxMult);
1772   fIntFlowCorrelationsAllVsMPro[5]->Sumw2();
1773   fIntFlowCorrelationsAllVsMPro[5]->GetXaxis()->SetTitle("M");  
1774   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[5]);
1775   fIntFlowCorrelationsAllVsMPro[6] = new TProfile("three3n2n1n","#LT#LT3#GT#GT_{3n|2n,n}",fnBinsMult,fMinMult,fMaxMult);
1776   fIntFlowCorrelationsAllVsMPro[6]->Sumw2();
1777   fIntFlowCorrelationsAllVsMPro[6]->GetXaxis()->SetTitle("M");  
1778   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[6]);
1779   fIntFlowCorrelationsAllVsMPro[7] = new TProfile("three4n2n2n","#LT#LT3#GT#GT_{4n|2n,2n}",fnBinsMult,fMinMult,fMaxMult);
1780   fIntFlowCorrelationsAllVsMPro[7]->Sumw2();
1781   fIntFlowCorrelationsAllVsMPro[7]->GetXaxis()->SetTitle("M");  
1782   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[7]);
1783   fIntFlowCorrelationsAllVsMPro[8] = new TProfile("three4n3n1n","#LT#LT3#GT#GT_{4n|3n,n}",fnBinsMult,fMinMult,fMaxMult);
1784   fIntFlowCorrelationsAllVsMPro[8]->Sumw2();
1785   fIntFlowCorrelationsAllVsMPro[8]->GetXaxis()->SetTitle("M");  
1786   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[8]);
1787   // 4-p correlations vs M:
1788   fIntFlowCorrelationsAllVsMPro[10] = new TProfile("four1n1n1n1n","#LT#LT4#GT#GT_{n,n|n,n}",fnBinsMult,fMinMult,fMaxMult);
1789   fIntFlowCorrelationsAllVsMPro[10]->Sumw2();
1790   fIntFlowCorrelationsAllVsMPro[10]->GetXaxis()->SetTitle("M");  
1791   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[10]);
1792   fIntFlowCorrelationsAllVsMPro[11] = new TProfile("four2n1n2n1n","#LT#LT4#GT#GT_{2n,n|2n,n}",fnBinsMult,fMinMult,fMaxMult);
1793   fIntFlowCorrelationsAllVsMPro[11]->Sumw2();
1794   fIntFlowCorrelationsAllVsMPro[11]->GetXaxis()->SetTitle("M");  
1795   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[11]);
1796   fIntFlowCorrelationsAllVsMPro[12] = new TProfile("four2n2n2n2n","#LT#LT4#GT#GT_{2n,2n|2n,2n}",fnBinsMult,fMinMult,fMaxMult);
1797   fIntFlowCorrelationsAllVsMPro[12]->Sumw2();
1798   fIntFlowCorrelationsAllVsMPro[12]->GetXaxis()->SetTitle("M");  
1799   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[12]);
1800   fIntFlowCorrelationsAllVsMPro[13] = new TProfile("four3n1n1n1n","#LT#LT4#GT#GT_{3n|n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1801   fIntFlowCorrelationsAllVsMPro[13]->Sumw2();
1802   fIntFlowCorrelationsAllVsMPro[13]->GetXaxis()->SetTitle("M");  
1803   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[13]);
1804   fIntFlowCorrelationsAllVsMPro[14] = new TProfile("four3n1n3n1n","#LT#LT4#GT#GT_{3n,n|3n,n}",fnBinsMult,fMinMult,fMaxMult);
1805   fIntFlowCorrelationsAllVsMPro[14]->Sumw2();
1806   fIntFlowCorrelationsAllVsMPro[14]->GetXaxis()->SetTitle("M");  
1807   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[14]);
1808   fIntFlowCorrelationsAllVsMPro[15] = new TProfile("four3n1n2n2n","#LT#LT4#GT#GT_{3n,n|2n,2n}",fnBinsMult,fMinMult,fMaxMult);
1809   fIntFlowCorrelationsAllVsMPro[15]->Sumw2();
1810   fIntFlowCorrelationsAllVsMPro[15]->GetXaxis()->SetTitle("M");  
1811   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[15]);
1812   fIntFlowCorrelationsAllVsMPro[16] = new TProfile("four4n2n1n1n","#LT#LT4#GT#GT_{4n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1813   fIntFlowCorrelationsAllVsMPro[16]->Sumw2();
1814   fIntFlowCorrelationsAllVsMPro[16]->GetXaxis()->SetTitle("M");  
1815   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[16]);
1816   // 5-p correlations vs M:
1817   fIntFlowCorrelationsAllVsMPro[18] = new TProfile("five2n1n1n1n1n","#LT#LT5#GT#GT_{2n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1818   fIntFlowCorrelationsAllVsMPro[18]->Sumw2();
1819   fIntFlowCorrelationsAllVsMPro[18]->GetXaxis()->SetTitle("M");  
1820   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[18]);
1821   fIntFlowCorrelationsAllVsMPro[19] = new TProfile("five2n2n2n1n1n","#LT#LT5#GT#GT_{2n,2n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1822   fIntFlowCorrelationsAllVsMPro[19]->Sumw2();
1823   fIntFlowCorrelationsAllVsMPro[19]->GetXaxis()->SetTitle("M");  
1824   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[19]);
1825   fIntFlowCorrelationsAllVsMPro[20] = new TProfile("five3n1n2n1n1n","#LT#LT5#GT#GT_{3n,n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1826   fIntFlowCorrelationsAllVsMPro[20]->Sumw2();
1827   fIntFlowCorrelationsAllVsMPro[20]->GetXaxis()->SetTitle("M");  
1828   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[20]);
1829   fIntFlowCorrelationsAllVsMPro[21] = new TProfile("five4n1n1n1n1n","#LT#LT5#GT#GT_{4n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1830   fIntFlowCorrelationsAllVsMPro[21]->Sumw2();
1831   fIntFlowCorrelationsAllVsMPro[21]->GetXaxis()->SetTitle("M");  
1832   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[21]);
1833   // 6-p correlations vs M:
1834   fIntFlowCorrelationsAllVsMPro[23] = new TProfile("six1n1n1n1n1n1n","#LT#LT6#GT#GT_{n,n,n|n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1835   fIntFlowCorrelationsAllVsMPro[23]->Sumw2();
1836   fIntFlowCorrelationsAllVsMPro[23]->GetXaxis()->SetTitle("M");  
1837   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[23]);
1838   fIntFlowCorrelationsAllVsMPro[24] = new TProfile("six2n1n1n2n1n1n","#LT#LT6#GT#GT_{2n,n,n|2n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1839   fIntFlowCorrelationsAllVsMPro[24]->Sumw2();
1840   fIntFlowCorrelationsAllVsMPro[24]->GetXaxis()->SetTitle("M");  
1841   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[24]);
1842   fIntFlowCorrelationsAllVsMPro[25] = new TProfile("six2n2n1n1n1n1n","#LT#LT6#GT#GT_{2n,2n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1843   fIntFlowCorrelationsAllVsMPro[25]->Sumw2();
1844   fIntFlowCorrelationsAllVsMPro[25]->GetXaxis()->SetTitle("M");  
1845   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[25]);
1846   fIntFlowCorrelationsAllVsMPro[26] = new TProfile("six3n1n1n1n1n1n","#LT#LT6#GT#GT_{3n,n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1847   fIntFlowCorrelationsAllVsMPro[26]->Sumw2();
1848   fIntFlowCorrelationsAllVsMPro[26]->GetXaxis()->SetTitle("M");  
1849   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[26]);
1850   // 7-p correlations vs M:
1851   fIntFlowCorrelationsAllVsMPro[28] = new TProfile("seven2n1n1n1n1n1n1n","#LT#LT7#GT#GT_{2n,n,n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1852   fIntFlowCorrelationsAllVsMPro[28]->Sumw2();
1853   fIntFlowCorrelationsAllVsMPro[28]->GetXaxis()->SetTitle("M");  
1854   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[28]);
1855   // 8-p correlations vs M:
1856   fIntFlowCorrelationsAllVsMPro[30] = new TProfile("eight1n1n1n1n1n1n1n1n","#LT#LT8#GT#GT_{n,n,n,n|n,n,n,n}",fnBinsMult,fMinMult,fMaxMult);
1857   fIntFlowCorrelationsAllVsMPro[30]->Sumw2();
1858   fIntFlowCorrelationsAllVsMPro[30]->GetXaxis()->SetTitle("M");  
1859   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[30]);
1860   // EXTRA correlations vs M for v3{5} study (to be improved - put them in a right order somewhere):
1861   fIntFlowCorrelationsAllVsMPro[32] = new TProfile("four4n2n3n3n","#LT#LT4#GT#GT_{4n,2n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
1862   fIntFlowCorrelationsAllVsMPro[32]->Sumw2();
1863   fIntFlowCorrelationsAllVsMPro[32]->GetXaxis()->SetTitle("M");  
1864   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[32]);
1865   fIntFlowCorrelationsAllVsMPro[33] = new TProfile("five3n3n2n2n2n","#LT#LT5#GT#GT_{3n,3n|2n,2n,2n}",fnBinsMult,fMinMult,fMaxMult);
1866   fIntFlowCorrelationsAllVsMPro[33]->Sumw2();
1867   fIntFlowCorrelationsAllVsMPro[33]->GetXaxis()->SetTitle("M");  
1868   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[33]);
1869   // EXTRA correlations vs M for Teaney-Yan study (to be improved - put them in a right order somewhere):
1870   fIntFlowCorrelationsAllVsMPro[34] = new TProfile("two5n5n","#LT#LT2#GT#GT_{5n|5n}",fnBinsMult,fMinMult,fMaxMult);
1871   fIntFlowCorrelationsAllVsMPro[34]->Sumw2();
1872   fIntFlowCorrelationsAllVsMPro[34]->GetXaxis()->SetTitle("M");  
1873   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[34]);  
1874   fIntFlowCorrelationsAllVsMPro[35] = new TProfile("two6n6n","#LT#LT2#GT#GT_{6n|6n}",fnBinsMult,fMinMult,fMaxMult);
1875   fIntFlowCorrelationsAllVsMPro[35]->Sumw2();
1876   fIntFlowCorrelationsAllVsMPro[35]->GetXaxis()->SetTitle("M");  
1877   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[35]);  
1878   fIntFlowCorrelationsAllVsMPro[36] = new TProfile("three5n3n2n","#LT#LT3#GT#GT_{5n|3n,2n}",fnBinsMult,fMinMult,fMaxMult);
1879   fIntFlowCorrelationsAllVsMPro[36]->Sumw2();
1880   fIntFlowCorrelationsAllVsMPro[36]->GetXaxis()->SetTitle("M");  
1881   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[36]); 
1882   fIntFlowCorrelationsAllVsMPro[37] = new TProfile("three5n4n1n","#LT#LT3#GT#GT_{5n|4n,1n}",fnBinsMult,fMinMult,fMaxMult);
1883   fIntFlowCorrelationsAllVsMPro[37]->Sumw2();
1884   fIntFlowCorrelationsAllVsMPro[37]->GetXaxis()->SetTitle("M");  
1885   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[37]);
1886   fIntFlowCorrelationsAllVsMPro[38] = new TProfile("three6n3n3n","#LT#LT3#GT#GT_{6n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
1887   fIntFlowCorrelationsAllVsMPro[38]->Sumw2();
1888   fIntFlowCorrelationsAllVsMPro[38]->GetXaxis()->SetTitle("M");  
1889   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[38]);  
1890   fIntFlowCorrelationsAllVsMPro[39] = new TProfile("three6n4n2n","#LT#LT3#GT#GT_{6n|4n,2n}",fnBinsMult,fMinMult,fMaxMult);
1891   fIntFlowCorrelationsAllVsMPro[39]->Sumw2();
1892   fIntFlowCorrelationsAllVsMPro[39]->GetXaxis()->SetTitle("M");  
1893   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[39]);
1894   fIntFlowCorrelationsAllVsMPro[40] = new TProfile("three6n5n1n","#LT#LT3#GT#GT_{6n|5n,1n}",fnBinsMult,fMinMult,fMaxMult);
1895   fIntFlowCorrelationsAllVsMPro[40]->Sumw2();
1896   fIntFlowCorrelationsAllVsMPro[40]->GetXaxis()->SetTitle("M");  
1897   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[40]);
1898   fIntFlowCorrelationsAllVsMPro[41] = new TProfile("four6n3n2n1n","#LT#LT4#GT#GT_{6n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
1899   fIntFlowCorrelationsAllVsMPro[41]->Sumw2();
1900   fIntFlowCorrelationsAllVsMPro[41]->GetXaxis()->SetTitle("M");  
1901   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[41]);
1902   fIntFlowCorrelationsAllVsMPro[42] = new TProfile("four3n2n3n2n","#LT#LT4#GT#GT_{3n,2n|3n,2n}",fnBinsMult,fMinMult,fMaxMult);
1903   fIntFlowCorrelationsAllVsMPro[42]->Sumw2();
1904   fIntFlowCorrelationsAllVsMPro[42]->GetXaxis()->SetTitle("M");  
1905   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[42]);
1906   fIntFlowCorrelationsAllVsMPro[43] = new TProfile("four4n1n3n2n","#LT#LT4#GT#GT_{4n,1n|3n,2n}",fnBinsMult,fMinMult,fMaxMult);
1907   fIntFlowCorrelationsAllVsMPro[43]->Sumw2();
1908   fIntFlowCorrelationsAllVsMPro[43]->GetXaxis()->SetTitle("M");  
1909   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[43]);
1910   fIntFlowCorrelationsAllVsMPro[44] = new TProfile("four3n3n3n3n","#LT#LT4#GT#GT_{3n,3n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
1911   fIntFlowCorrelationsAllVsMPro[44]->Sumw2();
1912   fIntFlowCorrelationsAllVsMPro[44]->GetXaxis()->SetTitle("M");  
1913   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[44]);
1914   fIntFlowCorrelationsAllVsMPro[45] = new TProfile("four4n2n3n3n","#LT#LT4#GT#GT_{4n,2n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
1915   fIntFlowCorrelationsAllVsMPro[45]->Sumw2();
1916   fIntFlowCorrelationsAllVsMPro[45]->GetXaxis()->SetTitle("M");  
1917   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[45]);
1918   fIntFlowCorrelationsAllVsMPro[46] = new TProfile("four5n1n3n3n","#LT#LT4#GT#GT_{5n,1n|3n,3n}",fnBinsMult,fMinMult,fMaxMult);
1919   fIntFlowCorrelationsAllVsMPro[46]->Sumw2();
1920   fIntFlowCorrelationsAllVsMPro[46]->GetXaxis()->SetTitle("M");  
1921   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[46]);
1922   fIntFlowCorrelationsAllVsMPro[47] = new TProfile("four4n2n4n2n","#LT#LT4#GT#GT_{4n,2n|4n,2n}",fnBinsMult,fMinMult,fMaxMult);
1923   fIntFlowCorrelationsAllVsMPro[47]->Sumw2();
1924   fIntFlowCorrelationsAllVsMPro[47]->GetXaxis()->SetTitle("M");  
1925   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[47]);
1926   fIntFlowCorrelationsAllVsMPro[48] = new TProfile("four5n1n4n2n","#LT#LT4#GT#GT_{5n,1n|4n,2n}",fnBinsMult,fMinMult,fMaxMult);
1927   fIntFlowCorrelationsAllVsMPro[48]->Sumw2();
1928   fIntFlowCorrelationsAllVsMPro[48]->GetXaxis()->SetTitle("M");  
1929   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[48]);
1930   fIntFlowCorrelationsAllVsMPro[49] = new TProfile("four5n3n1n1n","#LT#LT4#GT#GT_{5n|3n,1n,1n}",fnBinsMult,fMinMult,fMaxMult);
1931   fIntFlowCorrelationsAllVsMPro[49]->Sumw2();
1932   fIntFlowCorrelationsAllVsMPro[49]->GetXaxis()->SetTitle("M");  
1933   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[49]);
1934   fIntFlowCorrelationsAllVsMPro[50] = new TProfile("four5n2n2n1n","#LT#LT4#GT#GT_{5n|2n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
1935   fIntFlowCorrelationsAllVsMPro[50]->Sumw2();
1936   fIntFlowCorrelationsAllVsMPro[50]->GetXaxis()->SetTitle("M");  
1937   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[50]);
1938   fIntFlowCorrelationsAllVsMPro[51] = new TProfile("four5n1n5n1n","#LT#LT4#GT#GT_{5n,1n|5n,1n}",fnBinsMult,fMinMult,fMaxMult);
1939   fIntFlowCorrelationsAllVsMPro[51]->Sumw2();
1940   fIntFlowCorrelationsAllVsMPro[51]->GetXaxis()->SetTitle("M");  
1941   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[51]);
1942   fIntFlowCorrelationsAllVsMPro[52] = new TProfile("five3n3n3n2n1n","#LT#LT5#GT#GT_{3n,3n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
1943   fIntFlowCorrelationsAllVsMPro[52]->Sumw2();
1944   fIntFlowCorrelationsAllVsMPro[52]->GetXaxis()->SetTitle("M");  
1945   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[52]);
1946   fIntFlowCorrelationsAllVsMPro[53] = new TProfile("five4n2n3n2n1n","#LT#LT5#GT#GT_{4n,2n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
1947   fIntFlowCorrelationsAllVsMPro[53]->Sumw2();
1948   fIntFlowCorrelationsAllVsMPro[53]->GetXaxis()->SetTitle("M");  
1949   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[53]);
1950   fIntFlowCorrelationsAllVsMPro[54] = new TProfile("five3n2n3n1n1n","#LT#LT5#GT#GT_{3n,2n|3n,1n,1n}",fnBinsMult,fMinMult,fMaxMult);
1951   fIntFlowCorrelationsAllVsMPro[54]->Sumw2();
1952   fIntFlowCorrelationsAllVsMPro[54]->GetXaxis()->SetTitle("M");  
1953   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[54]);
1954   fIntFlowCorrelationsAllVsMPro[55] = new TProfile("five3n2n2n2n1n","#LT#LT5#GT#GT_{3n,2n|2n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
1955   fIntFlowCorrelationsAllVsMPro[55]->Sumw2();
1956   fIntFlowCorrelationsAllVsMPro[55]->GetXaxis()->SetTitle("M");  
1957   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[55]);
1958   fIntFlowCorrelationsAllVsMPro[56] = new TProfile("five5n1n3n2n1n","#LT#LT5#GT#GT_{5n,1n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
1959   fIntFlowCorrelationsAllVsMPro[56]->Sumw2();
1960   fIntFlowCorrelationsAllVsMPro[56]->GetXaxis()->SetTitle("M");  
1961   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[56]);
1962   fIntFlowCorrelationsAllVsMPro[57] = new TProfile("six3n2n1n3n2n1n","#LT#LT6#GT#GT_{3n,2n,1n|3n,2n,1n}",fnBinsMult,fMinMult,fMaxMult);
1963   fIntFlowCorrelationsAllVsMPro[57]->Sumw2();
1964   fIntFlowCorrelationsAllVsMPro[57]->GetXaxis()->SetTitle("M");  
1965   fIntFlowAllCorrelationsVsM->Add(fIntFlowCorrelationsAllVsMPro[57]);
1966  } // end of if(fCalculateAllCorrelationsVsM)
1967  // when particle weights are used some extra correlations appear:
1968  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights) 
1969  {
1970   TString intFlowExtraCorrelationsProName = "fIntFlowExtraCorrelationsPro";
1971   intFlowExtraCorrelationsProName += fAnalysisLabel->Data();
1972   fIntFlowExtraCorrelationsPro = new TProfile(intFlowExtraCorrelationsProName.Data(),"Average extra correlations for all events",100,0,100,"s");
1973   fIntFlowExtraCorrelationsPro->SetTickLength(-0.01,"Y");
1974   fIntFlowExtraCorrelationsPro->SetMarkerStyle(25);
1975   fIntFlowExtraCorrelationsPro->SetLabelSize(0.03);
1976   fIntFlowExtraCorrelationsPro->SetLabelOffset(0.01,"Y");
1977   // extra 2-p correlations:
1978   (fIntFlowExtraCorrelationsPro->GetXaxis())->SetBinLabel(1,"<<w1^3 w2 cos(n*(phi1-phi2))>>");
1979   (fIntFlowExtraCorrelationsPro->GetXaxis())->SetBinLabel(2,"<<w1 w2 w3^2 cos(n*(phi1-phi2))>>");
1980   fIntFlowProfiles->Add(fIntFlowExtraCorrelationsPro);
1981  } // end of if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
1982  // average product of correlations <2>, <4>, <6> and <8>:  
1983  TString productFlag[6] = {"<<2><4>>","<<2><6>>","<<2><8>>","<<4><6>>","<<4><8>>","<<6><8>>"};
1984  TString intFlowProductOfCorrelationsProName = "fIntFlowProductOfCorrelationsPro";
1985  intFlowProductOfCorrelationsProName += fAnalysisLabel->Data();
1986  fIntFlowProductOfCorrelationsPro = new TProfile(intFlowProductOfCorrelationsProName.Data(),"Average products of correlations",6,0,6);
1987  fIntFlowProductOfCorrelationsPro->SetTickLength(-0.01,"Y");
1988  fIntFlowProductOfCorrelationsPro->SetMarkerStyle(25); 
1989  fIntFlowProductOfCorrelationsPro->SetLabelSize(0.05);
1990  fIntFlowProductOfCorrelationsPro->SetLabelOffset(0.01,"Y");
1991  for(Int_t b=0;b<6;b++)
1992  {
1993   (fIntFlowProductOfCorrelationsPro->GetXaxis())->SetBinLabel(b+1,productFlag[b].Data());
1994  }
1995  fIntFlowProfiles->Add(fIntFlowProductOfCorrelationsPro); 
1996  // average product of correlations <2>, <4>, <6> and <8> versus multiplicity
1997  // [0=<<2><4>>,1=<<2><6>>,2=<<2><8>>,3=<<4><6>>,4=<<4><8>>,5=<<6><8>>]  
1998  if(fCalculateCumulantsVsM)
1999  {
2000   TString intFlowProductOfCorrelationsVsMProName = "fIntFlowProductOfCorrelationsVsMPro";
2001   intFlowProductOfCorrelationsVsMProName += fAnalysisLabel->Data();
2002   for(Int_t pi=0;pi<6;pi++)
2003   { 
2004    fIntFlowProductOfCorrelationsVsMPro[pi] = new TProfile(Form("%s, %s",intFlowProductOfCorrelationsVsMProName.Data(),productFlag[pi].Data()),
2005                                                           Form("%s versus multiplicity",productFlag[pi].Data()),
2006                                                           fnBinsMult,fMinMult,fMaxMult);             
2007    fIntFlowProductOfCorrelationsVsMPro[pi]->GetXaxis()->SetTitle("M");
2008    fIntFlowProfiles->Add(fIntFlowProductOfCorrelationsVsMPro[pi]);
2009   } // end of for(Int_t pi=0;pi<6;pi++)
2010  } // end of if(fCalculateCumulantsVsM) 
2011  // average product of correction terms for NUA:  
2012  TString intFlowProductOfCorrectionTermsForNUAProName = "fIntFlowProductOfCorrectionTermsForNUAPro";
2013  intFlowProductOfCorrectionTermsForNUAProName += fAnalysisLabel->Data();
2014  fIntFlowProductOfCorrectionTermsForNUAPro = new TProfile(intFlowProductOfCorrectionTermsForNUAProName.Data(),"Average products of correction terms for NUA",27,0,27);
2015  fIntFlowProductOfCorrectionTermsForNUAPro->SetTickLength(-0.01,"Y");
2016  fIntFlowProductOfCorrectionTermsForNUAPro->SetMarkerStyle(25); 
2017  fIntFlowProductOfCorrectionTermsForNUAPro->SetLabelSize(0.05);
2018  fIntFlowProductOfCorrectionTermsForNUAPro->SetLabelOffset(0.01,"Y");
2019  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(1,"<<2><cos(#phi)>>");
2020  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(2,"<<2><sin(#phi)>>");
2021  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(3,"<<cos(#phi)><sin(#phi)>>");
2022  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(4,"Cov(<2>,<cos(#phi_{1}+#phi_{2})>)");
2023  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(5,"Cov(<2>,<sin(#phi_{1}+#phi_{2})>)");
2024  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(6,"Cov(<2>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2025  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(7,"Cov(<2>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2026  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(8,"Cov(<4>,<cos(#phi)>)");
2027  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(9,"Cov(<4>,<sin(#phi)>)");
2028  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(10,"Cov(<4>,<cos(#phi_{1}+#phi_{2})>)");
2029  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(11,"Cov(<4>,<sin(#phi_{1}+#phi_{2})>)");
2030  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(12,"Cov(<4>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>>)");
2031  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(13,"Cov(<4>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>>)");
2032  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(14,"Cov(<cos(#phi)>,<cos(#phi_{1}+#phi_{2})>)"); 
2033  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(15,"Cov(<cos(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
2034  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(16,"Cov(<cos(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2035  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(17,"Cov(<cos(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2036  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(18,"Cov(<sin(#phi)>,<cos(#phi_{1}+#phi_{2})>)");
2037  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(19,"Cov(<sin(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
2038  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(20,"Cov(<sin(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2039  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(21,"Cov(<sin(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2040  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(22,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}+#phi_{2})>)");
2041  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(23,"Cov(<cos(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2042  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(24,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2043  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(25,"Cov(<sin(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2044  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(26,"Cov(<sin(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2045  (fIntFlowProductOfCorrectionTermsForNUAPro->GetXaxis())->SetBinLabel(27,"Cov(<cos(#phi_{1}-#phi_{2}-#phi_{3}>,<sin(#phi_{1}-#phi_{2}-#phi_{3}>)");
2046  fIntFlowProfiles->Add(fIntFlowProductOfCorrectionTermsForNUAPro);
2047  // average correction terms for non-uniform acceptance (with wrong errors!):
2048  for(Int_t sc=0;sc<2;sc++) // sin or cos terms
2049  {
2050   TString intFlowCorrectionTermsForNUAProName = "fIntFlowCorrectionTermsForNUAPro";
2051   intFlowCorrectionTermsForNUAProName += fAnalysisLabel->Data();
2052   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");
2053   fIntFlowCorrectionTermsForNUAPro[sc]->SetTickLength(-0.01,"Y");
2054   fIntFlowCorrectionTermsForNUAPro[sc]->SetMarkerStyle(25);
2055   fIntFlowCorrectionTermsForNUAPro[sc]->SetLabelSize(0.03);
2056   fIntFlowCorrectionTermsForNUAPro[sc]->SetLabelOffset(0.01,"Y");
2057   (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(1,Form("#LT#LT%s(n(phi1))#GT#GT",sinCosFlag[sc].Data()));
2058   (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(2,Form("#LT#LT%s(n(phi1+phi2))#GT#GT",sinCosFlag[sc].Data()));  
2059   (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(3,Form("#LT#LT%s(n(phi1-phi2-phi3))#GT#GT",sinCosFlag[sc].Data()));  
2060   (fIntFlowCorrectionTermsForNUAPro[sc]->GetXaxis())->SetBinLabel(4,Form("#LT#LT%s(n(2phi1-phi2))#GT#GT",sinCosFlag[sc].Data()));  
2061   fIntFlowProfiles->Add(fIntFlowCorrectionTermsForNUAPro[sc]);
2062   // versus multiplicity:
2063   if(fCalculateCumulantsVsM)
2064   {
2065    TString correctionTermFlag[4] = {"(n(phi1))","(n(phi1+phi2))","(n(phi1-phi2-phi3))","(n(2phi1-phi2))"}; // to be improved - hardwired 4
2066    for(Int_t ci=0;ci<4;ci++) // correction term index (to be improved - hardwired 4)
2067    {
2068     TString intFlowCorrectionTermsForNUAVsMProName = "fIntFlowCorrectionTermsForNUAVsMPro";
2069     intFlowCorrectionTermsForNUAVsMProName += fAnalysisLabel->Data();
2070     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");
2071     fIntFlowProfiles->Add(fIntFlowCorrectionTermsForNUAVsMPro[sc][ci]);
2072    }
2073   } // end of if(fCalculateCumulantsVsM)
2074  } // end of for(Int_t sc=0;sc<2;sc++) 
2075  
2076  // d) Book histograms holding the final results:
2077  // average correlations <<2>>, <<4>>, <<6>> and <<8>> for all events (with correct errors!):
2078  TString intFlowCorrelationsHistName = "fIntFlowCorrelationsHist";
2079  intFlowCorrelationsHistName += fAnalysisLabel->Data();
2080  fIntFlowCorrelationsHist = new TH1D(intFlowCorrelationsHistName.Data(),"Average correlations for all events",4,0,4);
2081  fIntFlowCorrelationsHist->SetTickLength(-0.01,"Y");
2082  fIntFlowCorrelationsHist->SetMarkerStyle(25);
2083  fIntFlowCorrelationsHist->SetLabelSize(0.06);
2084  fIntFlowCorrelationsHist->SetLabelOffset(0.01,"Y");
2085  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(1,"<<2>>");
2086  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(2,"<<4>>");
2087  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(3,"<<6>>");
2088  (fIntFlowCorrelationsHist->GetXaxis())->SetBinLabel(4,"<<8>>");
2089  fIntFlowResults->Add(fIntFlowCorrelationsHist);
2090  // average correlations <<2>>, <<4>>, <<6>> and <<8>> for all events (with correct errors!) vs M:
2091  if(fCalculateCumulantsVsM)
2092  {
2093   for(Int_t ci=0;ci<4;ci++) // correlation index
2094   {
2095    TString intFlowCorrelationsVsMHistName = "fIntFlowCorrelationsVsMHist";
2096    intFlowCorrelationsVsMHistName += fAnalysisLabel->Data();
2097    fIntFlowCorrelationsVsMHist[ci] = new TH1D(Form("%s, %s",intFlowCorrelationsVsMHistName.Data(),correlationFlag[ci].Data()),
2098                                               Form("%s vs multiplicity",correlationFlag[ci].Data()),
2099                                               fnBinsMult,fMinMult,fMaxMult);                                            
2100    fIntFlowCorrelationsVsMHist[ci]->GetYaxis()->SetTitle(correlationFlag[ci].Data());
2101    fIntFlowCorrelationsVsMHist[ci]->GetXaxis()->SetTitle("M");
2102    fIntFlowResults->Add(fIntFlowCorrelationsVsMHist[ci]);
2103   } // end of for(Int_t ci=0;ci<4;ci++) // correlation index   
2104  } // end of if(fCalculateCumulantsVsM) 
2105  // average all correlations for all events (with correct errors!):
2106  TString intFlowCorrelationsAllHistName = "fIntFlowCorrelationsAllHist";
2107  intFlowCorrelationsAllHistName += fAnalysisLabel->Data();
2108  fIntFlowCorrelationsAllHist = new TH1D(intFlowCorrelationsAllHistName.Data(),"Average correlations for all events",34,0,34);
2109  fIntFlowCorrelationsAllHist->SetTickLength(-0.01,"Y");
2110  fIntFlowCorrelationsAllHist->SetMarkerStyle(25);
2111  fIntFlowCorrelationsAllHist->SetLabelSize(0.03);
2112  fIntFlowCorrelationsAllHist->SetLabelOffset(0.01,"Y");
2113  // 2-p correlations:
2114  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
2115  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
2116  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
2117  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
2118  // 3-p correlations:
2119  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
2120  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
2121  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
2122  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
2123  // 4-p correlations:
2124  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
2125  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
2126  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
2127  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
2128  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
2129  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
2130  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
2131  // 5-p correlations:
2132  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
2133  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
2134  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
2135  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
2136  // 6-p correlations:
2137  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
2138  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
2139  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
2140  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
2141  // 7-p correlations:  
2142  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
2143  // 8-p correlations:
2144  (fIntFlowCorrelationsAllHist->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
2145  fIntFlowResults->Add(fIntFlowCorrelationsAllHist);
2146  // average correction terms for non-uniform acceptance (with correct errors!):
2147  for(Int_t sc=0;sc<2;sc++) // sin or cos terms
2148  {
2149   TString intFlowCorrectionTermsForNUAHistName = "fIntFlowCorrectionTermsForNUAHist";
2150   intFlowCorrectionTermsForNUAHistName += fAnalysisLabel->Data();
2151   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);
2152   fIntFlowCorrectionTermsForNUAHist[sc]->SetTickLength(-0.01,"Y");
2153   fIntFlowCorrectionTermsForNUAHist[sc]->SetMarkerStyle(25);
2154   fIntFlowCorrectionTermsForNUAHist[sc]->SetLabelSize(0.03);
2155   fIntFlowCorrectionTermsForNUAHist[sc]->SetLabelOffset(0.01,"Y");
2156   (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(1,Form("#LT#LT%s(n(#phi_{1}))#GT#GT",sinCosFlag[sc].Data()));
2157   (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(2,Form("#LT#LT%s(n(phi1+phi2))#GT#GT",sinCosFlag[sc].Data()));  
2158   (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(3,Form("#LT#LT%s(n(phi1-phi2-phi3))#GT#GT",sinCosFlag[sc].Data()));  
2159   (fIntFlowCorrectionTermsForNUAHist[sc]->GetXaxis())->SetBinLabel(4,Form("#LT#LT%s(n(2phi1-phi2))#GT#GT",sinCosFlag[sc].Data()));   
2160   fIntFlowResults->Add(fIntFlowCorrectionTermsForNUAHist[sc]);
2161  } // end of for(Int_t sc=0;sc<2;sc++) 
2162  // covariances (multiplied with weight dependent prefactor):
2163  TString intFlowCovariancesName = "fIntFlowCovariances";
2164  intFlowCovariancesName += fAnalysisLabel->Data();
2165  fIntFlowCovariances = new TH1D(intFlowCovariancesName.Data(),"Covariances (multiplied with weight dependent prefactor)",6,0,6);
2166  fIntFlowCovariances->SetLabelSize(0.04);
2167  fIntFlowCovariances->SetMarkerStyle(25);
2168  (fIntFlowCovariances->GetXaxis())->SetBinLabel(1,"Cov(<2>,<4>)");
2169  (fIntFlowCovariances->GetXaxis())->SetBinLabel(2,"Cov(<2>,<6>)");
2170  (fIntFlowCovariances->GetXaxis())->SetBinLabel(3,"Cov(<2>,<8>)");
2171  (fIntFlowCovariances->GetXaxis())->SetBinLabel(4,"Cov(<4>,<6>)");
2172  (fIntFlowCovariances->GetXaxis())->SetBinLabel(5,"Cov(<4>,<8>)");
2173  (fIntFlowCovariances->GetXaxis())->SetBinLabel(6,"Cov(<6>,<8>)");  
2174  fIntFlowResults->Add(fIntFlowCovariances);
2175  // sum of linear and quadratic event weights for <2>, <4>, <6> and <8>:
2176  TString intFlowSumOfEventWeightsName = "fIntFlowSumOfEventWeights";
2177  intFlowSumOfEventWeightsName += fAnalysisLabel->Data();
2178  for(Int_t power=0;power<2;power++)
2179  {
2180   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);
2181   fIntFlowSumOfEventWeights[power]->SetLabelSize(0.05);
2182   fIntFlowSumOfEventWeights[power]->SetMarkerStyle(25);
2183   if(power == 0)
2184   {
2185    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{<2>}");
2186    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{<4>}");
2187    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{<6>}");
2188    (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{<8>}");
2189   } else if (power == 1) 
2190     {
2191      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{<2>}^{2}");
2192      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{<4>}^{2}");
2193      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{<6>}^{2}");
2194      (fIntFlowSumOfEventWeights[power]->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{<8>}^{2}");
2195     }
2196   fIntFlowResults->Add(fIntFlowSumOfEventWeights[power]);
2197  } 
2198  // sum of products of event weights for correlations <2>, <4>, <6> and <8>:  
2199  TString intFlowSumOfProductOfEventWeightsName = "fIntFlowSumOfProductOfEventWeights";
2200  intFlowSumOfProductOfEventWeightsName += fAnalysisLabel->Data();
2201  fIntFlowSumOfProductOfEventWeights = new TH1D(intFlowSumOfProductOfEventWeightsName.Data(),"Sum of product of event weights for correlations",6,0,6);
2202  fIntFlowSumOfProductOfEventWeights->SetLabelSize(0.05);
2203  fIntFlowSumOfProductOfEventWeights->SetMarkerStyle(25);
2204  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{<2>} w_{<4>}");
2205  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{<2>} w_{<6>}");
2206  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{<2>} w_{<8>}");
2207  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(4,"#sum_{i=1}^{N} w_{<4>} w_{<6>}");
2208  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(5,"#sum_{i=1}^{N} w_{<4>} w_{<8>}");
2209  (fIntFlowSumOfProductOfEventWeights->GetXaxis())->SetBinLabel(6,"#sum_{i=1}^{N} w_{<6>} w_{<8>}");
2210  fIntFlowResults->Add(fIntFlowSumOfProductOfEventWeights);
2211  // final result for covariances of correlations (multiplied with weight dependent prefactor) versus M
2212  // [0=Cov(2,4),1=Cov(2,6),2=Cov(2,8),3=Cov(4,6),4=Cov(4,8),5=Cov(6,8)]:
2213  if(fCalculateCumulantsVsM)
2214  {
2215   TString intFlowCovariancesVsMName = "fIntFlowCovariancesVsM";
2216   intFlowCovariancesVsMName += fAnalysisLabel->Data();
2217   TString covarianceFlag[6] = {"Cov(<2>,<4>)","Cov(<2>,<6>)","Cov(<2>,<8>)","Cov(<4>,<6>)","Cov(<4>,<8>)","Cov(<6>,<8>)"};
2218   for(Int_t ci=0;ci<6;ci++)
2219   {
2220    fIntFlowCovariancesVsM[ci] = new TH1D(Form("%s, %s",intFlowCovariancesVsMName.Data(),covarianceFlag[ci].Data()),
2221                                          Form("%s vs multiplicity",covarianceFlag[ci].Data()),
2222                                          fnBinsMult,fMinMult,fMaxMult);
2223    fIntFlowCovariancesVsM[ci]->GetYaxis()->SetTitle(covarianceFlag[ci].Data());
2224    fIntFlowCovariancesVsM[ci]->GetXaxis()->SetTitle("M");
2225    fIntFlowResults->Add(fIntFlowCovariancesVsM[ci]);
2226   }
2227  } // end of if(fCalculateCumulantsVsM) 
2228  // sum of linear and quadratic event weights for <2>, <4>, <6> and <8> versus multiplicity
2229  // [0=sum{w_{<2>}},1=sum{w_{<4>}},2=sum{w_{<6>}},3=sum{w_{<8>}}][0=linear 1,1=quadratic]:
2230  if(fCalculateCumulantsVsM)
2231  {
2232   TString intFlowSumOfEventWeightsVsMName = "fIntFlowSumOfEventWeightsVsM";
2233   intFlowSumOfEventWeightsVsMName += fAnalysisLabel->Data();
2234   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>}"},
2235                            {"#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}"}};
2236   for(Int_t si=0;si<4;si++)
2237   {
2238    for(Int_t power=0;power<2;power++)
2239    {
2240     fIntFlowSumOfEventWeightsVsM[si][power] = new TH1D(Form("%s, %s",intFlowSumOfEventWeightsVsMName.Data(),sumFlag[power][si].Data()),
2241                                                        Form("%s vs multiplicity",sumFlag[power][si].Data()),
2242                                                        fnBinsMult,fMinMult,fMaxMult);    
2243     fIntFlowSumOfEventWeightsVsM[si][power]->GetYaxis()->SetTitle(sumFlag[power][si].Data());  
2244     fIntFlowSumOfEventWeightsVsM[si][power]->GetXaxis()->SetTitle("M");  
2245     fIntFlowResults->Add(fIntFlowSumOfEventWeightsVsM[si][power]);
2246    } // end of for(Int_t power=0;power<2;power++)
2247   } // end of for(Int_t si=0;si<4;si++)   
2248  } // end of if(fCalculateCumulantsVsM)
2249  // sum of products of event weights for correlations <2>, <4>, <6> and <8> vs M
2250  // [0=sum{w_{<2>}w_{<4>}},1=sum{w_{<2>}w_{<6>}},2=sum{w_{<2>}w_{<8>}},
2251  //  3=sum{w_{<4>}w_{<6>}},4=sum{w_{<4>}w_{<8>}},5=sum{w_{<6>}w_{<8>}}]:  
2252  if(fCalculateCumulantsVsM)
2253  {
2254   TString intFlowSumOfProductOfEventWeightsVsMName = "fIntFlowSumOfProductOfEventWeightsVsM";
2255   intFlowSumOfProductOfEventWeightsVsMName += fAnalysisLabel->Data();
2256   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>}",
2257                           "#sum_{i=1}^{N} w_{<4>} w_{<6>}","#sum_{i=1}^{N} w_{<4>} w_{<8>}","#sum_{i=1}^{N} w_{<6>} w_{<8>}"}; 
2258   for(Int_t pi=0;pi<6;pi++)
2259   {
2260    fIntFlowSumOfProductOfEventWeightsVsM[pi] = new TH1D(Form("%s, %s",intFlowSumOfProductOfEventWeightsVsMName.Data(),sopowFlag[pi].Data()),
2261                                                         Form("%s versus multiplicity",sopowFlag[pi].Data()),
2262                                                         fnBinsMult,fMinMult,fMaxMult); 
2263    fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetXaxis()->SetTitle("M");
2264    fIntFlowSumOfProductOfEventWeightsVsM[pi]->GetYaxis()->SetTitle(sopowFlag[pi].Data()); 
2265    fIntFlowResults->Add(fIntFlowSumOfProductOfEventWeightsVsM[pi]);
2266   } // end of for(Int_t pi=0;pi<6;pi++) 
2267  } // end of if(fCalculateCumulantsVsM)
2268  // covariances of NUA terms (multiplied with weight dependent prefactor):
2269  TString intFlowCovariancesNUAName = "fIntFlowCovariancesNUA";
2270  intFlowCovariancesNUAName += fAnalysisLabel->Data();
2271  fIntFlowCovariancesNUA = new TH1D(intFlowCovariancesNUAName.Data(),"Covariances for NUA (multiplied with weight dependent prefactor)",27,0,27);
2272  fIntFlowCovariancesNUA->SetLabelSize(0.04);
2273  fIntFlowCovariancesNUA->SetMarkerStyle(25);
2274  fIntFlowCovariancesNUA->GetXaxis()->SetLabelSize(0.02);
2275  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(1,"Cov(<2>,<cos(#phi)>");
2276  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(2,"Cov(<2>,<sin(#phi)>)");
2277  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(3,"Cov(<cos(#phi)>,<sin(#phi)>)");
2278  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(4,"Cov(<2>,<cos(#phi_{1}+#phi_{2})>)");
2279  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(5,"Cov(<2>,<sin(#phi_{1}+#phi_{2})>)");
2280  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(6,"Cov(<2>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2281  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(7,"Cov(<2>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2282  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(8,"Cov(<4>,<cos(#phi)>)");
2283  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(9,"Cov(<4>,<sin(#phi)>)");
2284  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(10,"Cov(<4>,<cos(#phi_{1}+#phi_{2})>)");
2285  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(11,"Cov(<4>,<sin(#phi_{1}+#phi_{2})>)");
2286  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(12,"Cov(<4>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>>)");
2287  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(13,"Cov(<4>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>>)");
2288  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(14,"Cov(<cos(#phi)>,<cos(#phi_{1}+#phi_{2})>)"); 
2289  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(15,"Cov(<cos(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
2290  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(16,"Cov(<cos(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2291  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(17,"Cov(<cos(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2292  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(18,"Cov(<sin(#phi)>,<cos(#phi_{1}+#phi_{2})>)");
2293  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(19,"Cov(<sin(#phi)>,<sin(#phi_{1}+#phi_{2})>)");
2294  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(20,"Cov(<sin(#phi)>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2295  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(21,"Cov(<sin(#phi)>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2296  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(22,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}+#phi_{2})>)");
2297  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(23,"Cov(<cos(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2298  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(24,"Cov(<cos(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2299  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(25,"Cov(<sin(#phi_{1}+#phi_{2})>,<cos(#phi_{1}-#phi_{2}-#phi_{3})>)");
2300  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(26,"Cov(<sin(#phi_{1}+#phi_{2})>,<sin(#phi_{1}-#phi_{2}-#phi_{3})>)");
2301  (fIntFlowCovariancesNUA->GetXaxis())->SetBinLabel(27,"Cov(<cos(#phi_{1}-#phi_{2}-#phi_{3}>,<sin(#phi_{1}-#phi_{2}-#phi_{3}>)");
2302  fIntFlowResults->Add(fIntFlowCovariancesNUA);
2303  // sum of linear and quadratic event weights for NUA terms:
2304  TString intFlowSumOfEventWeightsNUAName = "fIntFlowSumOfEventWeightsNUA";
2305  intFlowSumOfEventWeightsNUAName += fAnalysisLabel->Data();
2306  for(Int_t sc=0;sc<2;sc++)
2307  {
2308   for(Int_t power=0;power<2;power++)
2309   {
2310    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
2311    fIntFlowSumOfEventWeightsNUA[sc][power]->SetLabelSize(0.05);
2312    fIntFlowSumOfEventWeightsNUA[sc][power]->SetMarkerStyle(25);
2313    if(power == 0)
2314    {
2315     (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(1,Form("#sum_{i=1}^{N} w_{<%s(#phi)>}",sinCosFlag[sc].Data()));
2316     (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(2,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}+#phi_{2})>}",sinCosFlag[sc].Data()));
2317     (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(3,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}-#phi_{2}-#phi_{3})>}",sinCosFlag[sc].Data()));   
2318     (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(4,Form("#sum_{i=1}^{N} w_{<%s(2#phi_{1}-#phi_{2})>}",sinCosFlag[sc].Data()));
2319    } else if(power == 1) 
2320      {
2321       (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(1,Form("#sum_{i=1}^{N} w_{<%s(#phi)>}^{2}",sinCosFlag[sc].Data()));
2322       (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(2,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}+#phi_{2})>}^{2}",sinCosFlag[sc].Data()));
2323       (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(3,Form("#sum_{i=1}^{N} w_{<%s(#phi_{1}-#phi_{2}-#phi_{3})>}^{2}",sinCosFlag[sc].Data()));
2324       (fIntFlowSumOfEventWeightsNUA[sc][power]->GetXaxis())->SetBinLabel(4,Form("#sum_{i=1}^{N} w_{<%s(2#phi_{1}-#phi_{2})>}^{2}",sinCosFlag[sc].Data()));
2325      }
2326    fIntFlowResults->Add(fIntFlowSumOfEventWeightsNUA[sc][power]);
2327   }
2328  }  
2329  // sum of products of event weights for NUA terms:  
2330  TString intFlowSumOfProductOfEventWeightsNUAName = "fIntFlowSumOfProductOfEventWeightsNUA";
2331  intFlowSumOfProductOfEventWeightsNUAName += fAnalysisLabel->Data();
2332  fIntFlowSumOfProductOfEventWeightsNUA = new TH1D(intFlowSumOfProductOfEventWeightsNUAName.Data(),"Sum of product of event weights for NUA terms",27,0,27);
2333  fIntFlowSumOfProductOfEventWeightsNUA->SetLabelSize(0.05);
2334  fIntFlowSumOfProductOfEventWeightsNUA->SetMarkerStyle(25);
2335  (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(1,"#sum_{i=1}^{N} w_{<2>} w_{<cos(#phi)>}");
2336  (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(2,"#sum_{i=1}^{N} w_{<2>} w_{<sin(#phi)>}");
2337  (fIntFlowSumOfProductOfEventWeightsNUA->GetXaxis())->SetBinLabel(3,"#sum_{i=1}^{N} w_{<cos(#phi)>} w_{<sin(#phi)>}");
2338  // ....
2339  // to be improved - add labels for remaining bins
2340  // ....
2341  fIntFlowResults->Add(fIntFlowSumOfProductOfEventWeightsNUA);
2342  // Final results for reference Q-cumulants:
2343  TString cumulantFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"};
2344  TString intFlowQcumulantsName = "fIntFlowQcumulants";
2345  intFlowQcumulantsName += fAnalysisLabel->Data();
2346  fIntFlowQcumulants = new TH1D(intFlowQcumulantsName.Data(),"Reference Q-cumulants",4,0,4);
2347  if(fPropagateErrorAlsoFromNIT)
2348  {
2349   fIntFlowQcumulants->SetTitle("Reference Q-cumulants (error from non-isotropic terms also propagated)");
2350  }
2351  fIntFlowQcumulants->SetLabelSize(0.05);
2352  fIntFlowQcumulants->SetMarkerStyle(25);
2353  for(Int_t b=0;b<4;b++)
2354  {
2355   (fIntFlowQcumulants->GetXaxis())->SetBinLabel(b+1,cumulantFlag[b].Data());
2356  } 
2357  fIntFlowResults->Add(fIntFlowQcumulants);
2358  // Final results for reference Q-cumulants rebinned in M: 
2359  if(fCalculateCumulantsVsM)
2360  {
2361   TString intFlowQcumulantsRebinnedInMName = "fIntFlowQcumulantsRebinnedInM";
2362   intFlowQcumulantsRebinnedInMName += fAnalysisLabel->Data();
2363   fIntFlowQcumulantsRebinnedInM = new TH1D(intFlowQcumulantsRebinnedInMName.Data(),"Reference Q-cumulants rebinned in M",4,0,4);
2364   fIntFlowQcumulantsRebinnedInM->SetLabelSize(0.05);
2365   fIntFlowQcumulantsRebinnedInM->SetMarkerStyle(25);
2366   for(Int_t b=0;b<4;b++)
2367   {
2368    (fIntFlowQcumulantsRebinnedInM->GetXaxis())->SetBinLabel(b+1,cumulantFlag[b].Data());
2369   } 
2370   fIntFlowResults->Add(fIntFlowQcumulantsRebinnedInM);
2371  } // end of if(fCalculateCumulantsVsM) 
2372  // Ratio between error squared: with/without non-isotropic terms:
2373  TString intFlowQcumulantsErrorSquaredRatioName = "fIntFlowQcumulantsErrorSquaredRatio";
2374  intFlowQcumulantsErrorSquaredRatioName += fAnalysisLabel->Data();
2375  fIntFlowQcumulantsErrorSquaredRatio = new TH1D(intFlowQcumulantsErrorSquaredRatioName.Data(),"Error squared of reference Q-cumulants: #frac{with NUA terms}{without NUA terms}",4,0,4);
2376  fIntFlowQcumulantsErrorSquaredRatio->SetLabelSize(0.05);
2377  fIntFlowQcumulantsErrorSquaredRatio->SetMarkerStyle(25);
2378  for(Int_t b=0;b<4;b++)
2379  {
2380   (fIntFlowQcumulantsErrorSquaredRatio->GetXaxis())->SetBinLabel(b+1,cumulantFlag[b].Data());
2381  } 
2382  fIntFlowResults->Add(fIntFlowQcumulantsErrorSquaredRatio);
2383  // final results for integrated Q-cumulants versus multiplicity:
2384  if(fCalculateCumulantsVsM)
2385  {
2386   TString intFlowQcumulantsVsMName = "fIntFlowQcumulantsVsM";
2387   intFlowQcumulantsVsMName += fAnalysisLabel->Data();
2388   for(Int_t co=0;co<4;co++) // cumulant order
2389   {
2390    fIntFlowQcumulantsVsM[co] = new TH1D(Form("%s, %s",intFlowQcumulantsVsMName.Data(),cumulantFlag[co].Data()),
2391                                         Form("%s vs multipicity",cumulantFlag[co].Data()),
2392                                         fnBinsMult,fMinMult,fMaxMult);
2393    fIntFlowQcumulantsVsM[co]->GetXaxis()->SetTitle("M");                                     
2394    fIntFlowQcumulantsVsM[co]->GetYaxis()->SetTitle(cumulantFlag[co].Data());  
2395    fIntFlowResults->Add(fIntFlowQcumulantsVsM[co]);                                    
2396   } // end of for(Int_t co=0;co<4;co++) // cumulant order
2397  } // end of if(fCalculateCumulantsVsM)
2398  // final integrated flow estimates from Q-cumulants:
2399  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)};
2400  TString intFlowName = "fIntFlow";
2401  intFlowName += fAnalysisLabel->Data();  
2402  // integrated flow from Q-cumulants:
2403  fIntFlow = new TH1D(intFlowName.Data(),"Reference flow estimates from Q-cumulants",4,0,4);
2404  fIntFlow->SetLabelSize(0.05);
2405  fIntFlow->SetMarkerStyle(25);
2406  for(Int_t b=0;b<4;b++)
2407  {
2408   (fIntFlow->GetXaxis())->SetBinLabel(b+1,flowFlag[b].Data()); 
2409  }
2410  fIntFlowResults->Add(fIntFlow); 
2411  // Reference flow vs M rebinned in one huge bin:
2412  if(fCalculateCumulantsVsM)
2413  { 
2414   TString intFlowRebinnedInMName = "fIntFlowRebinnedInM";
2415   intFlowRebinnedInMName += fAnalysisLabel->Data();  
2416   fIntFlowRebinnedInM = new TH1D(intFlowRebinnedInMName.Data(),"Reference flow estimates from Q-cumulants (rebinned in M)",4,0,4);
2417   fIntFlowRebinnedInM->SetLabelSize(0.05);
2418   fIntFlowRebinnedInM->SetMarkerStyle(25);
2419   for(Int_t b=0;b<4;b++)
2420   {
2421    (fIntFlowRebinnedInM->GetXaxis())->SetBinLabel(b+1,flowFlag[b].Data()); 
2422   }
2423   fIntFlowResults->Add(fIntFlowRebinnedInM); 
2424  } 
2425  // integrated flow from Q-cumulants: versus multiplicity:
2426  if(fCalculateCumulantsVsM)
2427  {
2428   TString intFlowVsMName = "fIntFlowVsM";
2429   intFlowVsMName += fAnalysisLabel->Data();
2430   for(Int_t co=0;co<4;co++) // cumulant order
2431   {
2432    fIntFlowVsM[co] = new TH1D(Form("%s, %s",intFlowVsMName.Data(),flowFlag[co].Data()),
2433                               Form("%s vs multipicity",flowFlag[co].Data()),
2434                               fnBinsMult,fMinMult,fMaxMult);
2435    fIntFlowVsM[co]->GetXaxis()->SetTitle("M");                                     
2436    fIntFlowVsM[co]->GetYaxis()->SetTitle(flowFlag[co].Data());  
2437    fIntFlowResults->Add(fIntFlowVsM[co]);                                    
2438   } // end of for(Int_t co=0;co<4;co++) // cumulant order
2439  } // end of if(fCalculateCumulantsVsM)
2440  // quantifying detector effects effects to correlations:
2441  TString intFlowDetectorBiasName = "fIntFlowDetectorBias";
2442  intFlowDetectorBiasName += fAnalysisLabel->Data();  
2443  fIntFlowDetectorBias = new TH1D(intFlowDetectorBiasName.Data(),"Quantifying detector bias",4,0,4);
2444  fIntFlowDetectorBias->SetLabelSize(0.05);
2445  fIntFlowDetectorBias->SetMarkerStyle(25);
2446  for(Int_t ci=0;ci<4;ci++)
2447  {  
2448   (fIntFlowDetectorBias->GetXaxis())->SetBinLabel(ci+1,Form("#frac{corrected}{measured} %s",cumulantFlag[ci].Data()));
2449  }
2450  fIntFlowResults->Add(fIntFlowDetectorBias); 
2451  // quantifying detector effects to correlations versus multiplicity:
2452  if(fCalculateCumulantsVsM)
2453  {
2454   TString intFlowDetectorBiasVsMName = "fIntFlowDetectorBiasVsM";
2455   intFlowDetectorBiasVsMName += fAnalysisLabel->Data();
2456   for(Int_t ci=0;ci<4;ci++) // correlation index
2457   {
2458    fIntFlowDetectorBiasVsM[ci] = new TH1D(Form("%s for %s",intFlowDetectorBiasVsMName.Data(),cumulantFlag[ci].Data()),
2459                                           Form("Quantifying detector bias for %s vs multipicity",cumulantFlag[ci].Data()),
2460                                           fnBinsMult,fMinMult,fMaxMult);
2461    fIntFlowDetectorBiasVsM[ci]->GetXaxis()->SetTitle("M");                                     
2462    fIntFlowDetectorBiasVsM[ci]->GetYaxis()->SetTitle("#frac{corrected}{measured}");  
2463    fIntFlowResults->Add(fIntFlowDetectorBiasVsM[ci]);                                    
2464   } // end of for(Int_t co=0;co<4;co++) // cumulant order
2465  } // end of if(fCalculateCumulantsVsM)
2466    
2467 } // end of AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
2468
2469 //================================================================================================================================
2470
2471 void AliFlowAnalysisWithQCumulants::InitializeArraysForNestedLoops()
2472 {
2473  // Initialize arrays of all objects relevant for calculations with nested loops.
2474  
2475  // integrated flow:
2476  for(Int_t sc=0;sc<2;sc++) // sin or cos terms
2477  {
2478   fIntFlowDirectCorrectionTermsForNUA[sc] = NULL;
2479  } 
2480
2481  // differential flow:  
2482  // correlations:
2483  for(Int_t t=0;t<2;t++) // type: RP or POI
2484  { 
2485   for(Int_t pe=0;pe<2;pe++) // pt or eta
2486   {
2487    for(Int_t ci=0;ci<4;ci++) // correlation index
2488    {
2489     fDiffFlowDirectCorrelations[t][pe][ci] = NULL;
2490    } // end of for(Int_t ci=0;ci<4;ci++) // correlation index  
2491   } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
2492  } // end of for(Int_t t=0;t<2;t++) // type: RP or POI
2493  // correction terms for non-uniform acceptance:
2494  for(Int_t t=0;t<2;t++) // type: RP or 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 terms
2499    {
2500     for(Int_t cti=0;cti<9;cti++) // correction term index
2501     {
2502      fDiffFlowDirectCorrectionTermsForNUA[t][pe][sc][cti] = NULL;
2503     }   
2504    }
2505   } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
2506  } // end of for(Int_t t=0;t<2;t++) // type: RP or POI
2507
2508  // other differential correlators: 
2509  for(Int_t t=0;t<2;t++) // type: RP or POI
2510  { 
2511   for(Int_t pe=0;pe<2;pe++) // pt or eta
2512   {
2513    for(Int_t sc=0;sc<2;sc++) // sin or cos terms
2514    {
2515     for(Int_t ci=0;ci<1;ci++) // correlator index
2516     {
2517      fOtherDirectDiffCorrelators[t][pe][sc][ci] = NULL;
2518     }   
2519    }
2520   } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta
2521  } // end of for(Int_t t=0;t<2;t++) // type: RP or POI
2522
2523 } // end of void AliFlowAnalysisWithQCumulants::InitializeArraysForNestedLoops()
2524
2525 //================================================================================================================================
2526
2527 void AliFlowAnalysisWithQCumulants::BookEverythingForNestedLoops()
2528 {
2529  // Book all objects relevant for calculations with nested loops.
2530  
2531  TString sinCosFlag[2] = {"sin","cos"}; // to be improved (should I promote this to data members?)
2532  TString typeFlag[2] = {"RP","POI"}; // to be improved (should I promote this to data members?)
2533  TString ptEtaFlag[2] = {"p_{T}","#eta"}; // to be improved (should I promote this to data members?)
2534  TString reducedCorrelationIndex[4] = {"<2'>","<4'>","<6'>","<8'>"}; // to be improved (should I promote this to data members?)
2535  Double_t lowerPtEtaEdge[2] = {fPtMin+(fCrossCheckInPtBinNo-1)*fPtBinWidth,fEtaMin+(fCrossCheckInEtaBinNo-1)*fEtaBinWidth};
2536  Double_t upperPtEtaEdge[2] = {fPtMin+fCrossCheckInPtBinNo*fPtBinWidth,fEtaMin+fCrossCheckInEtaBinNo*fEtaBinWidth};
2537
2538  TString evaluateNestedLoopsName = "fEvaluateNestedLoops";
2539  evaluateNestedLoopsName += fAnalysisLabel->Data();
2540  fEvaluateNestedLoops = new TProfile(evaluateNestedLoopsName.Data(),"Flags for nested loops",4,0,4);
2541  fEvaluateNestedLoops->SetLabelSize(0.03);
2542  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(1,"fEvaluateIntFlowNestedLoops");
2543  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(2,"fEvaluateDiffFlowNestedLoops");
2544  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(3,"fCrossCheckInPtBinNo");
2545  (fEvaluateNestedLoops->GetXaxis())->SetBinLabel(4,"fCrossCheckInEtaBinNo");
2546  fEvaluateNestedLoops->Fill(0.5,(Int_t)fEvaluateIntFlowNestedLoops);
2547  fEvaluateNestedLoops->Fill(1.5,(Int_t)fEvaluateDiffFlowNestedLoops);
2548  fEvaluateNestedLoops->Fill(2.5,fCrossCheckInPtBinNo);
2549  fEvaluateNestedLoops->Fill(3.5,fCrossCheckInEtaBinNo);
2550  fNestedLoopsList->Add(fEvaluateNestedLoops);
2551  // nested loops for integrated flow:
2552  if(fEvaluateIntFlowNestedLoops)
2553  {
2554   // correlations:
2555   TString intFlowDirectCorrelationsName = "fIntFlowDirectCorrelations";
2556   intFlowDirectCorrelationsName += fAnalysisLabel->Data();
2557   fIntFlowDirectCorrelations = new TProfile(intFlowDirectCorrelationsName.Data(),"Multiparticle correlations calculated with nested loops (for int. flow)",58,0,58,"s");
2558   fNestedLoopsList->Add(fIntFlowDirectCorrelations);
2559   if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
2560   {
2561    TString intFlowExtraDirectCorrelationsName = "fIntFlowExtraDirectCorrelations";
2562    intFlowExtraDirectCorrelationsName += fAnalysisLabel->Data();
2563    fIntFlowExtraDirectCorrelations = new TProfile(intFlowExtraDirectCorrelationsName.Data(),"Extra multiparticle correlations calculated with nested loops (for int. flow)",100,0,100,"s");
2564    fNestedLoopsList->Add(fIntFlowExtraDirectCorrelations);  
2565   } // end of if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
2566   // correction terms for non-uniform acceptance:
2567   for(Int_t sc=0;sc<2;sc++) // sin or cos terms
2568   {
2569    TString intFlowDirectCorrectionTermsForNUAName = "fIntFlowDirectCorrectionTermsForNUA";
2570    intFlowDirectCorrectionTermsForNUAName += fAnalysisLabel->Data();
2571    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");
2572    fNestedLoopsList->Add(fIntFlowDirectCorrectionTermsForNUA[sc]);
2573   } // end of for(Int_t sc=0;sc<2;sc++) 
2574  } // end of if(fEvaluateIntFlowNestedLoops)
2575  
2576  // nested loops for differential flow: 
2577  if(fEvaluateDiffFlowNestedLoops)
2578  {
2579   // reduced correlations:
2580   TString diffFlowDirectCorrelationsName = "fDiffFlowDirectCorrelations";
2581   diffFlowDirectCorrelationsName += fAnalysisLabel->Data();
2582   for(Int_t t=0;t<2;t++) // type: RP or POI
2583   { 
2584    for(Int_t pe=0;pe<2;pe++) // pt or eta
2585    {
2586     for(Int_t rci=0;rci<4;rci++) // reduced correlation index
2587     {
2588      // reduced correlations:
2589      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");
2590      fDiffFlowDirectCorrelations[t][pe][rci]->SetXTitle(ptEtaFlag[pe].Data());
2591      fNestedLoopsList->Add(fDiffFlowDirectCorrelations[t][pe][rci]); // to be improved (add dedicated list to hold reduced correlations)
2592     } // end of for(Int_t rci=0;rci<4;rci++) // correlation index
2593    } // end of for(Int_t pe=0;pe<2;pe++) // pt or eta 
2594   } // end of for(Int_t t=0;t<2;t++) // type: RP or POI 
2595   
2596   
2597   // correction terms for non-uniform acceptance:
2598   TString diffFlowDirectCorrectionTermsForNUAName = "fDiffFlowDirectCorrectionTermsForNUA";
2599   diffFlowDirectCorrectionTermsForNUAName += fAnalysisLabel->Data();
2600   for(Int_t t=0;t<2;t++) // typeFlag (0 = RP, 1 = POI)
2601   { 
2602    for(Int_t pe=0;pe<2;pe++) // pt or eta
2603    {
2604     for(Int_t sc=0;sc<2;sc++) // sin or cos
2605     {
2606      for(Int_t cti=0;cti<9;cti++) // correction term index
2607      {
2608       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"); 
2609       fNestedLoopsList->Add(fDiffFlowDirectCorrectionTermsForNUA[t][pe][sc][cti]);
2610      }
2611     }
2612    }
2613   }
2614   // other differential correlators: 
2615   TString otherDirectDiffCorrelatorsName = "fOtherDirectDiffCorrelators";
2616   otherDirectDiffCorrelatorsName += fAnalysisLabel->Data();
2617   for(Int_t t=0;t<2;t++) // typeFlag (0 = RP, 1 = POI)
2618   { 
2619    for(Int_t pe=0;pe<2;pe++) // pt or eta
2620    {
2621     for(Int_t sc=0;sc<2;sc++) // sin or cos
2622     {
2623      for(Int_t ci=0;ci<1;ci++) // correlator index
2624      {
2625       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]); 
2626       fNestedLoopsList->Add(fOtherDirectDiffCorrelators[t][pe][sc][ci]);
2627      }
2628     }
2629    }
2630   }
2631   // number of RPs and POIs in selected pt and eta bins for cross-checkings:
2632   TString noOfParticlesInBinName = "fNoOfParticlesInBin";
2633   fNoOfParticlesInBin = new TH1D(noOfParticlesInBinName.Data(),"Number of RPs and POIs in selected p_{T} and #eta bin",4,0,4);
2634   fNoOfParticlesInBin->GetXaxis()->SetBinLabel(1,"# of RPs in p_{T} bin");
2635   fNoOfParticlesInBin->GetXaxis()->SetBinLabel(2,"# of RPs in #eta bin");
2636   fNoOfParticlesInBin->GetXaxis()->SetBinLabel(3,"# of POIs in p_{T} bin");
2637   fNoOfParticlesInBin->GetXaxis()->SetBinLabel(4,"# of POIs in #eta bin");
2638   fNestedLoopsList->Add(fNoOfParticlesInBin);
2639  } // end of if(fEvaluateDiffFlowNestedLoops)
2640
2641 } // end of AliFlowAnalysisWithQCumulants::BookEverythingForNestedLoops()
2642
2643 //================================================================================================================================
2644
2645 void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations()
2646 {
2647  // Calculate in this method all multiparticle azimuthal correlations.
2648  //
2649  // Remark 1: All multiparticle correlations are stored in TProfile fIntFlowCorrelationsAllPro;
2650  // Remark 2: There is a special TProfile fIntFlowCorrelationsPro holding results 
2651  //           only for same harmonic's correlations <<2>>, <<4>>, <<6>> and <<8>>;  
2652  // Remark 3: Binning of fIntFlowCorrelationsAllPro is organized as follows:
2653  // --------------------------------------------------------------------------------------------------------------------
2654  //  1st bin: <2>_{1n|1n} = two1n1n = cos(1n(phi1-phi2))>
2655  //  2nd bin: <2>_{2n|2n} = two2n2n = cos(2n(phi1-phi2))>
2656  //  3rd bin: <2>_{3n|3n} = two3n3n = cos(3n(phi1-phi2))> 
2657  //  4th bin: <2>_{4n|4n} = two4n4n = cos(4n(phi1-phi2))>
2658  //  5th bin:           ----  EMPTY ----
2659  //  6th bin: <3>_{2n|1n,1n} = three2n1n1n = <cos(n(2*phi1-phi2-phi3))>
2660  //  7th bin: <3>_{3n|2n,1n} = three3n2n1n = <cos(n(3*phi1-2*phi2-phi3))>
2661  //  8th bin: <3>_{4n|2n,2n} = three4n2n2n = <cos(n(4*phi1-2*phi2-2*phi3))>
2662  //  9th bin: <3>_{4n|3n,1n} = three4n3n1n = <cos(n(4*phi1-3*phi2-phi3))>
2663  // 10th bin:           ----  EMPTY ----
2664  // 11th bin: <4>_{1n,1n|1n,1n} = four1n1n1n1n = <cos(n(phi1+phi2-phi3-phi4))>
2665  // 12th bin: <4>_{2n,1n|2n,1n} = four2n1n2n1n = <cos(n(2*phi1+phi2-2*phi3-phi4))>
2666  // 13th bin: <4>_{2n,2n|2n,2n} = four2n2n2n2n = <cos(2n(phi1+phi2-phi3-phi4))>
2667  // 14th bin: <4>_{3n|1n,1n,1n} = four3n1n1n1n = <cos(n(3*phi1-phi2-phi3-phi4))> 
2668  // 15th bin: <4>_{3n,1n|3n,1n} = four3n1n3n1n = <cos(n(3*phi1+phi2-3*phi3-phi4))>
2669  // 16th bin: <4>_{3n,1n|2n,2n} = four3n1n2n2n = <cos(n(3*phi1+phi2-2*phi3-2*phi4))>
2670  // 17th bin: <4>_{4n|2n,1n,1n} = four4n2n1n1n = <cos(n(4*phi1-2*phi2-phi3-phi4))>
2671  // 18th bin:           ----  EMPTY ----
2672  // 19th bin: <5>_{2n,1n|1n,1n,1n} = five2n1n1n1n1n = <cos(n(2*phi1+phi2-phi3-phi4-phi5))>
2673  // 20th bin: <5>_{2n,2n|2n,1n,1n} = five2n2n2n1n1n = <cos(n(2*phi1+2*phi2-2*phi3-phi4-phi5))>
2674  // 21st bin: <5>_{3n,1n|2n,1n,1n} = five3n1n2n1n1n = <cos(n(3*phi1+phi2-2*phi3-phi4-phi5))>
2675  // 22nd bin: <5>_{4n|1n,1n,1n,1n} = five4n1n1n1n1n = <cos(n(4*phi1-phi2-phi3-phi4-phi5))>
2676  // 23rd bin:           ----  EMPTY ----
2677  // 24th bin: <6>_{1n,1n,1n|1n,1n,1n} = six1n1n1n1n1n1n = <cos(n(phi1+phi2+phi3-phi4-phi5-phi6))>
2678  // 25th bin: <6>_{2n,1n,1n|2n,1n,1n} = six2n1n1n2n1n1n = <cos(n(2*phi1+phi2+phi3-2*phi4-phi5-phi6))>
2679  // 26th bin: <6>_{2n,2n|1n,1n,1n,1n} = six2n2n1n1n1n1n = <cos(n(2*phi1+2*phi2-phi3-phi4-phi5-phi6))>
2680  // 27th bin: <6>_{3n,1n|1n,1n,1n,1n} = six3n1n1n1n1n1n = <cos(n(3*phi1+phi2-phi3-phi4-phi5-phi6))> 
2681  // 28th bin:           ----  EMPTY ----
2682  // 29th bin: <7>_{2n,1n,1n|1n,1n,1n,1n} = seven2n1n1n1n1n1n1n =  <cos(n(2*phi1+phi2+phi3-phi4-phi5-phi6-phi7))>
2683  // 30th bin:           ----  EMPTY ----
2684  // 31st bin: <8>_{1n,1n,1n,1n|1n,1n,1n,1n} = eight1n1n1n1n1n1n1n1n = <cos(n(phi1+phi2+phi3+phi4-phi5-phi6-phi7-phi8))>
2685  // 32nd bin:           ----  EMPTY ----
2686  //  Extra correlations for v3{5} study: 
2687  // 33rd bin: <4>_{4n,2n|3n,3n} = four4n2n3n3n = <cos(n(4*phi1+2*phi2-3*phi3-3*phi4))>
2688  // 34th bin: <5>_{3n,3n|2n,2n,2n} = five3n3n2n2n2n = <cos(n(3*phi1+3*phi2-2*phi3-2*phi4-2*phi5))> 
2689  //  Extra correlations for Teaney-Yan study: 
2690  // 35th bin: <2>_{5n|5n} = two5n5n = <cos(5n(phi1-phi2)> 
2691  // 36th bin: <2>_{6n|6n} = two6n6n = <cos(6n(phi1-phi2)> 
2692  // 37th bin: <3>_{5n|3n,2n} = three5n3n2n = <cos(n(5*phi1-3*phi2-2*phi3)> 
2693  // 38th bin: <3>_{5n|4n,1n} = three5n4n1n = <cos(n(5*phi1-4*phi2-1*phi3)> 
2694  // 39th bin: <3>_{6n|3n,3n} = three6n3n3n = <cos(n(6*phi1-3*phi2-3*phi3)> 
2695  // 40th bin: <3>_{6n|4n,2n} = three6n4n2n = <cos(n(6*phi1-4*phi2-2*phi3)> 
2696  // 41st bin: <3>_{6n|5n,1n} = three6n5n1n = <cos(n(6*phi1-5*phi2-1*phi3)>
2697  // 42nd bin: <4>_{6n|3n,2n,1n} = four6n3n2n1n = <cos(n(6*phi1-3*phi2-2*phi3-1*phi4)>
2698  // 43rd bin: <4>_{3n,2n|3n,2n} = four3n2n3n2n = <cos(n(3*phi1+2*phi2-3*phi3-2*phi4)>
2699  // 44th bin: <4>_{4n,1n|3n,2n} = four4n1n3n2n = <cos(n(4*phi1+1*phi2-3*phi3-2*phi4)>
2700  // 45th bin: <4>_{3n,3n|3n,3n} = four3n3n3n3n = <cos(3n*(phi1+phi2-phi3-phi4))> 
2701  // 46th bin: <4>_{4n,2n|3n,3n} = four4n2n3n3n = <cos(n(4*phi1+2*phi2-3*phi3-3*phi4)>
2702  // 47th bin: <4>_{5n,1n|3n,3n} = four5n1n3n3n = <cos(n(5*phi1+1*phi2-3*phi3-3*phi4)>
2703  // 48th bin: <4>_{4n,2n|4n,2n} = four4n2n4n2n = <cos(n(4*phi1+2*phi2-4*phi3-2*phi4)> 
2704  // 49th bin: <4>_{5n,1n|4n,2n} = four5n1n4n2n = <cos(n(5*phi1+1*phi2-4*phi3-2*phi4)>
2705  // 50th bin: <4>_{5n|3n,1n,1n} = four5n3n1n1n = <cos(n(5*phi1-3*phi2-1*phi3-1*phi4)>
2706  // 51st bin: <4>_{5n|2n,2n,1n} = four5n2n2n1n = <cos(n(5*phi1-2*phi2-2*phi3-1*phi4)>
2707  // 52nd bin: <4>_{5n,1n|5n,1n} = four5n1n5n1n = <cos(n(5*phi1+1*phi2-5*phi3-1*phi4)>
2708  // 53rd bin: <5>_{3n,3n|3n,2n,1n} = five3n3n3n2n1n = <cos(n(3*phi1+3*phi2-3*phi3-2*phi4-1*phi5)>
2709  // 54th bin: <5>_{4n,2n|3n,2n,1n} = five4n2n3n2n1n = <cos(n(4*phi1+2*phi2-3*phi3-2*phi4-1*phi5)>
2710  // 55th bin: <5>_{3n,2n|3n,1n,1n} = five3n2n3n1n1n = <cos(n(3*phi1+2*phi2-3*phi3-1*phi4-1*phi5)>
2711  // 56th bin: <5>_{3n,2n|2n,2n,1n} = five3n2n2n2n1n = <cos(n(3*phi1+2*phi2-2*phi3-2*phi4-1*phi5)>
2712  // 57th bin: <5>_{5n,1n|3n,2n,1n} = five5n1n3n2n1n = <cos(n(5*phi1+1*phi2-3*phi3-2*phi4-1*phi5)>
2713  // 58th bin: <6>_{3n,2n,1n|3n,2n,1n} = six3n2n1n3n2n1n = <cos(n(3*phi1+2*phi2+1*phi3-3*phi4-2*phi5-1*phi6)>
2714  // --------------------------------------------------------------------------------------------------------------------
2715  
2716  // Multiplicity of an event:
2717  Double_t dMult = (*fSpk)(0,0);
2718  // Real parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n, 4n, 5n and 6n: 
2719  Double_t dReQ1n = (*fReQ)(0,0);
2720  Double_t dReQ2n = (*fReQ)(1,0);
2721  Double_t dReQ3n = (*fReQ)(2,0);
2722  Double_t dReQ4n = (*fReQ)(3,0);
2723  Double_t dReQ5n = (*fReQ)(4,0); 
2724  Double_t dReQ6n = (*fReQ)(5,0);
2725  // Imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n, 4n, 5n and 6n:
2726  Double_t dImQ1n = (*fImQ)(0,0);
2727  Double_t dImQ2n = (*fImQ)(1,0);
2728  Double_t dImQ3n = (*fImQ)(2,0);
2729  Double_t dImQ4n = (*fImQ)(3,0);
2730  Double_t dImQ5n = (*fImQ)(4,0); 
2731  Double_t dImQ6n = (*fImQ)(5,0);
2732   
2733  // Real parts of expressions involving various combinations of Q-vectors which appears
2734  // simultaneously in several equations for multiparticle correlations bellow: 
2735  // Re[Q_{2n}Q_{n}^*Q_{n}^*]
2736  Double_t reQ2nQ1nstarQ1nstar = pow(dReQ1n,2.)*dReQ2n+2.*dReQ1n*dImQ1n*dImQ2n-pow(dImQ1n,2.)*dReQ2n; 
2737  // Re[Q_{6n}Q_{3n}^*Q_{3n}^*]
2738  Double_t reQ6nQ3nstarQ3nstar = pow(dReQ3n,2.)*dReQ6n+2.*dReQ3n*dImQ3n*dImQ6n-pow(dImQ3n,2.)*dReQ6n;  
2739  // Re[Q_{4n}Q_{2n}^*Q_{2n}^*]
2740  Double_t reQ4nQ2nstarQ2nstar = pow(dReQ2n,2.)*dReQ4n+2.*dReQ2n*dImQ2n*dImQ4n-pow(dImQ2n,2.)*dReQ4n;
2741  // Re[Q_{4n}Q_{3n}^*Q_{n}^*]
2742  Double_t reQ4nQ3nstarQ1nstar = dReQ4n*(dReQ3n*dReQ1n-dImQ3n*dImQ1n)+dImQ4n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);
2743  // Re[Q_{3n}Q_{2n}^*Q_{n}^*]
2744  Double_t reQ3nQ2nstarQ1nstar = dReQ3n*dReQ2n*dReQ1n-dReQ3n*dImQ2n*dImQ1n+dImQ3n*dReQ2n*dImQ1n
2745                               + dImQ3n*dImQ2n*dReQ1n; 
2746  // Re[Q_{5n}Q_{3n}^*Q_{2n}^*]
2747  Double_t reQ5nQ3nstarQ2nstar = dReQ5n*dReQ2n*dReQ3n-dReQ5n*dImQ2n*dImQ3n+dImQ5n*dReQ2n*dImQ3n
2748                               + dImQ5n*dImQ2n*dReQ3n;                             
2749  // Re[Q_{5n}Q_{4n}^*Q_{1n}^*]
2750  Double_t reQ5nQ4nstarQ1nstar = dReQ5n*dReQ4n*dReQ1n-dReQ5n*dImQ4n*dImQ1n+dImQ5n*dReQ4n*dImQ1n
2751                               + dImQ5n*dImQ4n*dReQ1n;                              
2752  // Re[Q_{6n}Q_{5n}^*Q_{1n}^*]                              
2753  Double_t reQ6nQ5nstarQ1nstar = dReQ6n*dReQ5n*dReQ1n-dReQ6n*dImQ5n*dImQ1n+dImQ6n*dReQ5n*dImQ1n
2754                               + dImQ6n*dImQ5n*dReQ1n;
2755  // Re[Q_{6n}Q_{4n}^*Q_{2n}^*]                              
2756  Double_t reQ6nQ4nstarQ2nstar = dReQ6n*dReQ4n*dReQ2n-dReQ6n*dImQ4n*dImQ2n+dImQ6n*dReQ4n*dImQ2n
2757                               + dImQ6n*dImQ4n*dReQ2n;
2758  // Re[Q_{3n}Q_{n}Q_{2n}^*Q_{2n}^*]
2759  Double_t reQ3nQ1nQ2nstarQ2nstar = (pow(dReQ2n,2.)-pow(dImQ2n,2.))*(dReQ3n*dReQ1n-dImQ3n*dImQ1n) 
2760                                  + 2.*dReQ2n*dImQ2n*(dReQ3n*dImQ1n+dImQ3n*dReQ1n);                                
2761  // Re[Q_{3n}Q_{n}^*Q_{n}^*Q_{n}^*]
2762  Double_t reQ3nQ1nstarQ1nstarQ1nstar = dReQ3n*pow(dReQ1n,3)-3.*dReQ1n*dReQ3n*pow(dImQ1n,2)
2763                                      + 3.*dImQ1n*dImQ3n*pow(dReQ1n,2)-dImQ3n*pow(dImQ1n,3);
2764  // Re[Q_{4n}Q_{2n}^*Q_{n}^*Q_{n}^*]
2765  Double_t reQ4nQ2nstarQ1nstarQ1nstar = (dReQ4n*dReQ2n+dImQ4n*dImQ2n)*(pow(dReQ1n,2)-pow(dImQ1n,2)) 
2766                                      + 2.*dReQ1n*dImQ1n*(dImQ4n*dReQ2n-dReQ4n*dImQ2n);  
2767  // Re[Q_{4n}Q_{2n}^*Q_{3n}^*Q_{3n}^*]
2768  Double_t reQ4nQ2nQ3nstarQ3nstar = (dReQ4n*dReQ2n-dImQ4n*dImQ2n)*(dReQ3n*dReQ3n-dImQ3n*dImQ3n)
2769                                  + 2.*(dReQ4n*dImQ2n+dImQ4n*dReQ2n)*dReQ3n*dImQ3n;                    
2770  // Re[Q_{4n}Q_{n}Q_{3n}^*Q_{2n}^*]
2771  Double_t reQ4nQ1nQ3nstarQ2nstar = dImQ1n*dImQ2n*dImQ3n*dImQ4n+dImQ3n*dImQ4n*dReQ1n*dReQ2n 
2772                                  + dImQ2n*dImQ4n*dReQ1n*dReQ3n-dImQ1n*dImQ4n*dReQ2n*dReQ3n
2773                                  - dImQ2n*dImQ3n*dReQ1n*dReQ4n+dImQ1n*dImQ3n*dReQ2n*dReQ4n 
2774                                  + dImQ1n*dImQ2n*dReQ3n*dReQ4n+dReQ1n*dReQ2n*dReQ3n*dReQ4n;
2775  // Re[Q_{5n}Q_{n}Q_{4n}^*Q_{2n}^*]
2776  Double_t reQ5nQ1nQ4nstarQ2nstar = dImQ1n*dImQ2n*dImQ4n*dImQ5n+dImQ4n*dImQ5n*dReQ1n*dReQ2n 
2777                                  + dImQ2n*dImQ5n*dReQ1n*dReQ4n-dImQ1n*dImQ5n*dReQ2n*dReQ4n
2778                                  - dImQ2n*dImQ4n*dReQ1n*dReQ5n+dImQ1n*dImQ4n*dReQ2n*dReQ5n 
2779                                  + dImQ1n*dImQ2n*dReQ4n*dReQ5n+dReQ1n*dReQ2n*dReQ4n*dReQ5n;                                  
2780  // Re[Q_{5n}Q_{n}Q_{3n}^*Q_{3n}^*]                                  
2781  Double_t reQ5nQ1nQ3nstarQ3nstar = dImQ1n*pow(dImQ3n,2.)*dImQ5n+2.*dImQ3n*dImQ5n*dReQ1n*dReQ3n
2782                                  - dImQ1n*dImQ5n*pow(dReQ3n,2.)-pow(dImQ3n,2.)*dReQ1n*dReQ5n 
2783                                  + 2.*dImQ1n*dImQ3n*dReQ3n*dReQ5n+dReQ1n*pow(dReQ3n,2.)*dReQ5n;
2784  // Re[Q_{5n}Q_{3n}^*Q_{n}^*Q_{n}^*]                                  
2785  Double_t reQ5nQ3nstarQ1nstarQ1nstar = -pow(dImQ1n,2.)*dImQ3n*dImQ5n+dImQ3n*dImQ5n*pow(dReQ1n,2.)
2786                                      + 2.*dImQ1n*dImQ5n*dReQ1n*dReQ3n-2.*dImQ1n*dImQ3n*dReQ1n*dReQ5n 
2787                                      - pow(dImQ1n,2.)*dReQ3n*dReQ5n+pow(dReQ1n,2.)*dReQ3n*dReQ5n;                     
2788  // Re[Q_{5n}Q_{2n}^*Q_{2n}^*Q_{n}^*]                                  
2789  Double_t reQ5nQ2nstarQ2nstarQ1nstar = -pow(dImQ2n,2.)*dImQ1n*dImQ5n+dImQ1n*dImQ5n*pow(dReQ2n,2.)
2790                                      + 2.*dImQ2n*dImQ5n*dReQ2n*dReQ1n-2.*dImQ2n*dImQ1n*dReQ2n*dReQ5n 
2791                                      - pow(dImQ2n,2.)*dReQ1n*dReQ5n+pow(dReQ2n,2.)*dReQ1n*dReQ5n;
2792  // Re[Q_{6n}Q_{2n}^*Q_{2n}^*Q_{2n}^*]
2793  Double_t reQ6nQ2nstarQ2nstarQ2nstar = dReQ6n*pow(dReQ2n,3)-3.*dReQ2n*dReQ6n*pow(dImQ2n,2)
2794                                      + 3.*dImQ2n*dImQ6n*pow(dReQ2n,2)-dImQ6n*pow(dImQ2n,3);                                        
2795  // |Q_{2n}|^2 |Q_{n}|^2
2796  Double_t dQ2nQ1nQ2nstarQ1nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))*(pow(dReQ1n,2.)+pow(dImQ1n,2.));
2797  // |Q_{4n}|^2 |Q_{2n}|^2
2798  Double_t dQ4nQ2nQ4nstarQ2nstar = (pow(dReQ4n,2.)+pow(dImQ4n,2.))*(pow(dReQ2n,2.)+pow(dImQ2n,2.));
2799  // |Q_{3n}|^2 |Q_{2n}|^2
2800  Double_t dQ3nQ2nQ3nstarQ2nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))*(pow(dReQ3n,2.)+pow(dImQ3n,2.));
2801  // |Q_{5n}|^2 |Q_{n}|^2
2802  Double_t dQ5nQ1nQ5nstarQ1nstar = (pow(dReQ5n,2.)+pow(dImQ5n,2.))*(pow(dReQ1n,2.)+pow(dImQ1n,2.));
2803  // Re[Q_{2n}Q_{n}Q_{n}^*Q_{n}^*Q_{n}^*]
2804  Double_t reQ2nQ1nQ1nstarQ1nstarQ1nstar = (dReQ2n*dReQ1n-dImQ2n*dImQ1n)*(pow(dReQ1n,3)-3.*dReQ1n*pow(dImQ1n,2))
2805                                         + (dReQ2n*dImQ1n+dReQ1n*dImQ2n)*(3.*dImQ1n*pow(dReQ1n,2)-pow(dImQ1n,3)); 
2806  // Re[Q_{2n}Q_{2n}Q_{2n}^*Q_{n}^*Q_{n}^*]
2807  Double_t reQ2nQ2nQ2nstarQ1nstarQ1nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))
2808                                         * (dReQ2n*(pow(dReQ1n,2.)-pow(dImQ1n,2.)) + 2.*dImQ2n*dReQ1n*dImQ1n);
2809  // Re[Q_{4n}Q_{n}^*Q_{n}^*Q_{n}^*Q_{n}^*]
2810  Double_t reQ4nQ1nstarQ1nstarQ1nstarQ1nstar = pow(dReQ1n,4.)*dReQ4n-6.*pow(dReQ1n,2.)*dReQ4n*pow(dImQ1n,2.)
2811                                             + pow(dImQ1n,4.)*dReQ4n+4.*pow(dReQ1n,3.)*dImQ1n*dImQ4n
2812                                             - 4.*pow(dImQ1n,3.)*dReQ1n*dImQ4n;
2813  // Re[Q_{3n}Q_{n}Q_{2n}^*Q_{n}^*Q_{n}^*]
2814  Double_t reQ3nQ1nQ2nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)+pow(dImQ1n,2.))
2815                                         * (dReQ1n*dReQ2n*dReQ3n-dReQ3n*dImQ1n*dImQ2n
2816                                         + dReQ2n*dImQ1n*dImQ3n+dReQ1n*dImQ2n*dImQ3n);
2817  // Re[Q_{6n}Q_{n}Q_{3n}^*Q_{2n}^*Q_{n}^*]
2818  Double_t reQ6nQ3nstarQ2nstarQ1nstar = dReQ1n*dReQ2n*dReQ3n*dReQ6n-dReQ3n*dReQ6n*dImQ1n*dImQ2n
2819                                      - dReQ2n*dReQ6n*dImQ1n*dImQ3n-dReQ1n*dReQ6n*dImQ2n*dImQ3n
2820                                      + dReQ2n*dReQ3n*dImQ1n*dImQ6n+dReQ1n*dReQ3n*dImQ2n*dImQ6n 
2821                                      + dReQ1n*dReQ2n*dImQ3n*dImQ6n-dImQ1n*dImQ2n*dImQ3n*dImQ6n;
2822  // Re[Q_{3n}Q_{3n}Q_{3n}^*Q_{2n}^*Q_{n}^*]
2823  Double_t reQ3nQ3nQ3nstarQ2nstarQ1nstar = (pow(dImQ3n,2.)+pow(dReQ3n,2.))
2824                                         * (dImQ2n*dImQ3n*dReQ1n+dImQ1n*dImQ3n*dReQ2n
2825                                         - dImQ1n*dImQ2n*dReQ3n+dReQ1n*dReQ2n*dReQ3n);   
2826  // Re[Q_{3n}Q_{3n}Q_{2n}^*Q_{2n}^*Q_{2n}^*]
2827  Double_t reQ3nQ3nQ2nstarQ2nstarQ2nstar = pow(dReQ2n,3.)*pow(dReQ3n,2.) 
2828                                         - 3.*dReQ2n*pow(dReQ3n,2.)*pow(dImQ2n,2.)
2829                                         + 6.*pow(dReQ2n,2.)*dReQ3n*dImQ2n*dImQ3n 
2830                                         - 2.*dReQ3n*pow(dImQ2n,3.)*dImQ3n-pow(dReQ2n,3.)*pow(dImQ3n,2.) 
2831                                         + 3.*dReQ2n*pow(dImQ2n,2.)*pow(dImQ3n,2.);
2832  // Re[Q_{4n}Q_{2n}Q_{3n}^*Q_{2n}^*Q_{n}^*]
2833  Double_t reQ4nQ2nQ3nstarQ2nstarQ1nstar = (pow(dImQ2n,2.)+pow(dReQ2n,2.))
2834                                         * (dImQ3n*dImQ4n*dReQ1n+dImQ1n*dImQ4n*dReQ3n 
2835                                         - dImQ1n*dImQ3n*dReQ4n+dReQ1n*dReQ3n*dReQ4n);
2836  // Re[Q_{3n}Q_{2n}Q_{3n}^*Q_{n}^*Q_{n}^*]
2837  Double_t reQ3nQ2nQ3nstarQ1nstarQ1nstar = -(pow(dImQ3n,2.)+pow(dReQ3n,2.))
2838                                         * (-2.*dImQ1n*dImQ2n*dReQ1n+pow(dImQ1n,2.)*dReQ2n-pow(dReQ1n,2.)*dReQ2n);                              
2839  // Re[Q_{3n}Q_{2n}Q_{2n}^*Q_{2n}^*Q_{n}^*]
2840  Double_t reQ3nQ2nQ2nstarQ2nstarQ1nstar = (pow(dImQ2n,2.)+pow(dReQ2n,2.))
2841                                         * (dImQ2n*dImQ3n*dReQ1n+dImQ1n*dImQ3n*dReQ2n 
2842                                         - dImQ1n*dImQ2n*dReQ3n+dReQ1n*dReQ2n*dReQ3n);
2843  // Re[Q_{5n}Q_{n}Q_{3n}^*Q_{2n}^*Q_{n}^*]
2844  Double_t reQ5nQ1nQ3nstarQ2nstarQ1nstar = (pow(dImQ1n,2.)+pow(dReQ1n,2.))
2845                                         * (dImQ3n*dImQ5n*dReQ2n+dImQ2n*dImQ5n*dReQ3n 
2846                                         - dImQ2n*dImQ3n*dReQ5n+dReQ2n*dReQ3n*dReQ5n);   
2847  // Re[Q_{2n}Q_{2n}Q_{n}^*Q_{n}^*Q_{n}^*Q_{n}^*]
2848  Double_t reQ2nQ2nQ1nstarQ1nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)*dReQ2n-2.*dReQ1n*dReQ2n*dImQ1n-dReQ2n*pow(dImQ1n,2.)
2849                                                + dImQ2n*pow(dReQ1n,2.)+2.*dReQ1n*dImQ1n*dImQ2n-pow(dImQ1n,2.)*dImQ2n)
2850                                                * (pow(dReQ1n,2.)*dReQ2n+2.*dReQ1n*dReQ2n*dImQ1n-dReQ2n*pow(dImQ1n,2.)
2851                                                - dImQ2n*pow(dReQ1n,2.)+2.*dReQ1n*dImQ1n*dImQ2n+pow(dImQ1n,2.)*dImQ2n); 
2852  // Re[Q_{3n}Q_{n}Q_{n}^*Q_{n}^*Q_{n}^*Q_{n}^*]
2853  Double_t reQ3nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = (pow(dReQ1n,2.)+pow(dImQ1n,2.))
2854                                                * (pow(dReQ1n,3.)*dReQ3n-3.*dReQ1n*dReQ3n*pow(dImQ1n,2.)
2855                                                + 3.*pow(dReQ1n,2.)*dImQ1n*dImQ3n-pow(dImQ1n,3.)*dImQ3n);
2856  // |Q_{2n}|^2 |Q_{n}|^4
2857  Double_t dQ2nQ1nQ1nQ2nstarQ1nstarQ1nstar = (pow(dReQ2n,2.)+pow(dImQ2n,2.))*pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.); 
2858  // |Q_{3n}|^2 |Q_{2n}|^2 |Q_{n}|^2
2859  Double_t dQ3nQ2nQ1nQ3nstarQ2nstarQ1nstar = (pow(dReQ3n,2.)+pow(dImQ3n,2.))*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
2860                                           * (pow(dReQ1n,2.)+pow(dImQ1n,2.));
2861  // Re[Q_{2n}Q_{n}Q_{n}Q_{n}^*Q_{n}^*Q_{n}^*Q_{n}^*]
2862  Double_t reQ2nQ1nQ1nQ1nstarQ1nstarQ1nstarQ1nstar = pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)
2863                                                   * (pow(dReQ1n,2.)*dReQ2n-dReQ2n*pow(dImQ1n,2.)
2864                                                   + 2.*dReQ1n*dImQ1n*dImQ2n);
2865     
2866  // Results for multiparticle azimuthal correlations:
2867  // 2-particle:
2868  Double_t two1n1n = 0.; // <cos(n(phi1-phi2))>
2869  Double_t two2n2n = 0.; // <cos(2n(phi1-phi2))>
2870  Double_t two3n3n = 0.; // <cos(3n(phi1-phi2))>
2871  Double_t two4n4n = 0.; // <cos(4n(phi1-phi2))>
2872  if(dMult>1)
2873  {
2874   two1n1n = (pow(dReQ1n,2.)+pow(dImQ1n,2.)-dMult)/(dMult*(dMult-1.)); 
2875   two2n2n = (pow(dReQ2n,2.)+pow(dImQ2n,2.)-dMult)/(dMult*(dMult-1.)); 
2876   two3n3n = (pow(dReQ3n,2.)+pow(dImQ3n,2.)-dMult)/(dMult*(dMult-1.)); 
2877   two4n4n = (pow(dReQ4n,2.)+pow(dImQ4n,2.)-dMult)/(dMult*(dMult-1.)); 
2878   // Average 2-particle correlations for single event: 
2879   fIntFlowCorrelationsAllEBE->SetBinContent(1,two1n1n);
2880   fIntFlowCorrelationsAllEBE->SetBinContent(2,two2n2n);
2881   fIntFlowCorrelationsAllEBE->SetBinContent(3,two3n3n);
2882   fIntFlowCorrelationsAllEBE->SetBinContent(4,two4n4n);         
2883   // Average 2-particle correlations for all events:      
2884   fIntFlowCorrelationsAllPro->Fill(0.5,two1n1n,dMult*(dMult-1.));
2885   fIntFlowCorrelationsAllPro->Fill(1.5,two2n2n,dMult*(dMult-1.)); 
2886   fIntFlowCorrelationsAllPro->Fill(2.5,two3n3n,dMult*(dMult-1.)); 
2887   fIntFlowCorrelationsAllPro->Fill(3.5,two4n4n,dMult*(dMult-1.)); 
2888   // Store separetately <2>:
2889   fIntFlowCorrelationsEBE->SetBinContent(1,two1n1n); // <2>  
2890   // Testing other multiplicity weights:
2891   Double_t mWeight2p = 0.;
2892   if(!strcmp(fMultiplicityWeight->Data(),"combinations"))
2893   {
2894    mWeight2p = dMult*(dMult-1.);
2895   } else if(!strcmp(fMultiplicityWeight->Data(),"unit"))
2896     {
2897      mWeight2p = 1.;    
2898     } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity"))
2899       {
2900        mWeight2p = dMult;           
2901       }          
2902   fIntFlowEventWeightsForCorrelationsEBE->SetBinContent(1,mWeight2p); // eW_<2>
2903   fIntFlowCorrelationsPro->Fill(0.5,two1n1n,mWeight2p);
2904   fIntFlowSquaredCorrelationsPro->Fill(0.5,two1n1n*two1n1n,mWeight2p);
2905   if(fCalculateCumulantsVsM)
2906   {
2907    fIntFlowCorrelationsVsMPro[0]->Fill(dMult+0.5,two1n1n,mWeight2p);
2908    fIntFlowSquaredCorrelationsVsMPro[0]->Fill(dMult+0.5,two1n1n*two1n1n,mWeight2p);
2909   } 
2910   if(fCalculateAllCorrelationsVsM)
2911   {
2912    fIntFlowCorrelationsAllVsMPro[0]->Fill(dMult+0.5,two1n1n,mWeight2p);
2913    fIntFlowCorrelationsAllVsMPro[1]->Fill(dMult+0.5,two2n2n,mWeight2p);
2914    fIntFlowCorrelationsAllVsMPro[2]->Fill(dMult+0.5,two3n3n,mWeight2p);
2915    fIntFlowCorrelationsAllVsMPro[3]->Fill(dMult+0.5,two4n4n,mWeight2p);
2916   }  
2917  } // end of if(dMult>1)
2918  
2919  // 3-particle:
2920  Double_t three2n1n1n = 0.; // <cos(n(2*phi1-phi2-phi3))>
2921  Double_t three3n2n1n = 0.; // <cos(n(3*phi1-2*phi2-phi3))>
2922  Double_t three4n2n2n = 0.; // <cos(n(4*phi1-2*phi2-2*phi3))>
2923  Double_t three4n3n1n = 0.; // <cos(n(4*phi1-3*phi2-phi3))> 
2924  if(dMult>2)
2925  {
2926   three2n1n1n = (reQ2nQ1nstarQ1nstar-2.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))
2927               - (pow(dReQ2n,2.)+pow(dImQ2n,2.))+2.*dMult)
2928               / (dMult*(dMult-1.)*(dMult-2.));                     
2929   three3n2n1n = (reQ3nQ2nstarQ1nstar-(pow(dReQ3n,2.)+pow(dImQ3n,2.))
2930               - (pow(dReQ2n,2.)+pow(dImQ2n,2.))
2931               - (pow(dReQ1n,2.)+pow(dImQ1n,2.))+2.*dMult)
2932               / (dMult*(dMult-1.)*(dMult-2.));
2933   three4n2n2n = (reQ4nQ2nstarQ2nstar-2.*(pow(dReQ2n,2.)+pow(dImQ2n,2.))
2934               - (pow(dReQ4n,2.)+pow(dImQ4n,2.))+2.*dMult)
2935               / (dMult*(dMult-1.)*(dMult-2.)); 
2936   three4n3n1n = (reQ4nQ3nstarQ1nstar-(pow(dReQ4n,2.)+pow(dImQ4n,2.))
2937               - (pow(dReQ3n,2.)+pow(dImQ3n,2.))
2938               - (pow(dReQ1n,2.)+pow(dImQ1n,2.))+2.*dMult)
2939               / (dMult*(dMult-1.)*(dMult-2.));              
2940   // Average 3-particle correlations for single event: 
2941   fIntFlowCorrelationsAllEBE->SetBinContent(6,three2n1n1n);
2942   fIntFlowCorrelationsAllEBE->SetBinContent(7,three3n2n1n);
2943   fIntFlowCorrelationsAllEBE->SetBinContent(8,three4n2n2n);
2944   fIntFlowCorrelationsAllEBE->SetBinContent(9,three4n3n1n);
2945   // Average 3-particle correlations for all events:                
2946   fIntFlowCorrelationsAllPro->Fill(5.5,three2n1n1n,dMult*(dMult-1.)*(dMult-2.)); 
2947   fIntFlowCorrelationsAllPro->Fill(6.5,three3n2n1n,dMult*(dMult-1.)*(dMult-2.));
2948   fIntFlowCorrelationsAllPro->Fill(7.5,three4n2n2n,dMult*(dMult-1.)*(dMult-2.)); 
2949   fIntFlowCorrelationsAllPro->Fill(8.5,three4n3n1n,dMult*(dMult-1.)*(dMult-2.));  
2950   // Average 3-particle correlations vs M for all events:                
2951   if(fCalculateAllCorrelationsVsM)
2952   {
2953    fIntFlowCorrelationsAllVsMPro[5]->Fill(dMult+0.5,three2n1n1n,dMult*(dMult-1.)*(dMult-2.));
2954    fIntFlowCorrelationsAllVsMPro[6]->Fill(dMult+0.5,three3n2n1n,dMult*(dMult-1.)*(dMult-2.));
2955    fIntFlowCorrelationsAllVsMPro[7]->Fill(dMult+0.5,three4n2n2n,dMult*(dMult-1.)*(dMult-2.));
2956    fIntFlowCorrelationsAllVsMPro[8]->Fill(dMult+0.5,three4n3n1n,dMult*(dMult-1.)*(dMult-2.));
2957   }    
2958  } // end of if(dMult>2)
2959  
2960  // 4-particle:
2961  Double_t four1n1n1n1n = 0.; // <cos(n(phi1+phi2-phi3-phi4))>
2962  Double_t four2n2n2n2n = 0.; // <cos(2n(phi1+phi2-phi3-phi4))>
2963  Double_t four2n1n2n1n = 0.; // <cos(n(2*phi1+phi2-2*phi3-phi4))> 
2964  Double_t four3n1n1n1n = 0.; // <cos(n(3*phi1-phi2-phi3-phi4))> 
2965  Double_t four4n2n1n1n = 0.; // <cos(n(4*phi1-2*phi2-phi3-phi4))> 
2966  Double_t four3n1n2n2n = 0.; // <cos(n(3*phi1+phi2-2*phi3-2*phi4))> 
2967  Double_t four3n1n3n1n = 0.; // <cos(n(3*phi1+phi2-3*phi3-phi4))>    
2968  if(dMult>3)
2969  {
2970   four1n1n1n1n = (2.*dMult*(dMult-3.)+pow((pow(dReQ1n,2.)+pow(dImQ1n,2.)),2.)-4.*(dMult-2.)*(pow(dReQ1n,2.)