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