Changes to compile with Root6
[u/mrichter/AliRoot.git] / PWG / FLOW / Base / AliFlowAnalysisWithMultiparticleCorrelations.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 multi-particle *
18  *           correlations            * 
19  *                                   * 
20  * author: Ante Bilandzic            * 
21  *        (abilandzic@gmail.com)     *
22  ************************************/ 
23
24 #define AliFlowAnalysisWithMultiparticleCorrelations_cxx
25
26 #include "AliFlowAnalysisWithMultiparticleCorrelations.h"
27
28 using std::endl;
29 using std::cout;
30 using std::flush;
31 using std::ofstream;
32 using std::ios;
33
34 //================================================================================================================
35
36 ClassImp(AliFlowAnalysisWithMultiparticleCorrelations)
37
38 AliFlowAnalysisWithMultiparticleCorrelations::AliFlowAnalysisWithMultiparticleCorrelations(): 
39  // 0.) Base:
40  fHistList(NULL),
41  fInternalFlagsPro(NULL),
42  fUseInternalFlags(kFALSE),
43  fMinNoRPs(-44),
44  fMaxNoRPs(-44),
45  fExactNoRPs(-44),
46  fPropagateError(kTRUE),
47  fAnalysisTag(""),
48  fDumpThePoints(kFALSE),
49  fMaxNoEventsPerFile(100),
50  // 1.) Control histograms:
51  fControlHistogramsList(NULL),
52  fControlHistogramsFlagsPro(NULL),
53  fFillControlHistograms(kFALSE),
54  fFillKinematicsHist(kFALSE),
55  fFillMultDistributionsHist(kFALSE),   
56  fFillMultCorrelationsHist(kFALSE),
57  // 2.) Q-vector:
58  fQvectorList(NULL),       
59  fQvectorFlagsPro(NULL),
60  fCalculateQvector(kFALSE),
61  fCalculateDiffQvectors(kFALSE),
62  // 3.) Correlations:
63  fCorrelationsList(NULL),
64  fCorrelationsFlagsPro(NULL),
65  fCalculateCorrelations(kFALSE),
66  fMaxHarmonic(6), // TBI this shall not be hardwired in the ideal world...
67  fMaxCorrelator(8),
68  fCalculateIsotropic(kFALSE),
69  fCalculateSame(kFALSE),
70  fSkipZeroHarmonics(kFALSE),
71  fCalculateSameIsotropic(kFALSE),
72  fCalculateAll(kFALSE),
73  fDontGoBeyond(0),
74  fCalculateOnlyForHarmonicQC(kFALSE),
75  fCalculateOnlyForSC(kFALSE),
76  fCalculateOnlyCos(kFALSE),
77  fCalculateOnlySin(kFALSE),
78  // 4.) Event-by-event cumulants:
79  fEbECumulantsList(NULL),
80  fEbECumulantsFlagsPro(NULL),
81  fCalculateEbECumulants(kFALSE),
82  // 5.) Weights:
83  fWeightsList(NULL),
84  fWeightsFlagsPro(NULL),
85  // 6.) Nested loops:
86  fNestedLoopsList(NULL),
87  fNestedLoopsFlagsPro(NULL),
88  fCrossCheckWithNestedLoops(kFALSE),
89  fCrossCheckDiffWithNestedLoops(kFALSE),
90  fNestedLoopsResultsCosPro(NULL),
91  fNestedLoopsResultsSinPro(NULL),
92  fNestedLoopsDiffResultsPro(NULL),
93  // 7.) 'Standard candles':
94  fStandardCandlesList(NULL),
95  fStandardCandlesFlagsPro(NULL),
96  fCalculateStandardCandles(kFALSE),
97  fPropagateErrorSC(kTRUE),
98  fStandardCandlesHist(NULL),
99  fProductsSCPro(NULL), 
100  // 8.) Q-cumulants:
101  fQcumulantsList(NULL),
102  fQcumulantsFlagsPro(NULL),
103  fCalculateQcumulants(kFALSE),
104  fHarmonicQC(2),
105  fPropagateErrorQC(kTRUE),
106  fQcumulantsHist(NULL),
107  fReferenceFlowHist(NULL),
108  fProductsQCPro(NULL),
109  // 9.) Differential correlations:
110  fDiffCorrelationsList(NULL),
111  fDiffCorrelationsFlagsPro(NULL),
112  fCalculateDiffCorrelations(kFALSE),
113  fCalculateDiffCos(kTRUE),
114  fCalculateDiffSin(kFALSE),
115  fCalculateDiffCorrelationsVsPt(kTRUE),
116  fUseDefaultBinning(kTRUE),
117  fnDiffBins(-44),
118  fRangesDiffBins(NULL),
119  fDiffBinNo(-1)
120  {
121   // Constructor.  
122   
123   // a) Book grandmother of all lists;
124   // b) Initialize all arrays.
125
126   // a) Book grandmother of all lists:
127   fHistList = new TList();
128   fHistList->SetName("cobjMPC");
129   fHistList->SetOwner(kTRUE);
130
131   // b) Initialize all arrays:
132   this->InitializeArraysForControlHistograms();
133   this->InitializeArraysForQvector();
134   this->InitializeArraysForCorrelations();
135   this->InitializeArraysForEbECumulants();
136   this->InitializeArraysForWeights();
137   this->InitializeArraysForQcumulants();
138   this->InitializeArraysForDiffCorrelations(); 
139   this->InitializeArraysForNestedLoops(); 
140
141  } // end of AliFlowAnalysisWithMultiparticleCorrelations::AliFlowAnalysisWithMultiparticleCorrelations()
142  
143 //================================================================================================================  
144
145 AliFlowAnalysisWithMultiparticleCorrelations::~AliFlowAnalysisWithMultiparticleCorrelations()
146 {
147  // Destructor.
148  
149  delete fHistList;
150
151 } // end of AliFlowAnalysisWithMultiparticleCorrelations::~AliFlowAnalysisWithMultiparticleCorrelations()
152
153 //================================================================================================================
154
155 void AliFlowAnalysisWithMultiparticleCorrelations::Init()
156 {
157  // Well, this is method Init().
158
159  // a) Trick to avoid name clashes, part 1; 
160  // b) Cross-check the initial settings before starting this adventure;
161  // c) Book all objects;
162  // d) Set all flags;
163  // *) Trick to avoid name clashes, part 2. 
164
165  // a) Trick to avoid name clashes, part 1: 
166  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus(); 
167  TH1::AddDirectory(kFALSE);
168  
169  // b) Cross-check the initial settings before starting this adventure:
170  this->CrossCheckSettings();
171
172  // c) Book all objects:
173  this->BookAndNestAllLists();
174  this->BookEverythingForBase();
175  this->BookEverythingForControlHistograms();
176  this->BookEverythingForQvector();
177  this->BookEverythingForWeights();
178  this->BookEverythingForCorrelations();
179  this->BookEverythingForEbECumulants();
180  this->BookEverythingForNestedLoops();
181  this->BookEverythingForStandardCandles();
182  this->BookEverythingForQcumulants();
183  this->BookEverythingForDiffCorrelations(); 
184
185  // d) Set all flags:
186  // ... 
187
188  // *) Trick to avoid name clashes, part 2:
189  TH1::AddDirectory(oldHistAddStatus);
190
191 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::Init()
192
193 //================================================================================================================
194
195 void AliFlowAnalysisWithMultiparticleCorrelations::Make(AliFlowEventSimple *anEvent)
196 {
197  // Running over data only in this method.
198  
199  // a) Cross-check internal flags;
200  // b) Cross-check all pointers used in this method;
201  // c) Fill control histograms;
202  // d) Fill Q-vector components;
203  // e) Calculate multi-particle correlations from Q-vector components; 
204  // f) Calculate e-b-e cumulants; 
205  // g) Reset Q-vector components;
206  // h) Cross-check results with nested loops;
207  // i) Dump the points.
208
209  // a) Cross-check internal flags:
210  if(fUseInternalFlags){if(!this->CrossCheckInternalFlags(anEvent)){return;}}
211
212  // b) Cross-check all pointers used in this method:
213  this->CrossCheckPointersUsedInMake(); // TBI shall I call this method first  
214
215  // c) Fill control histograms:
216  if(fFillControlHistograms){this->FillControlHistograms(anEvent);}
217  
218  // d) Fill Q-vector components:
219  if(fCalculateQvector||fCalculateDiffQvectors){this->FillQvector(anEvent);}
220
221  // e) Calculate multi-particle correlations from Q-vector components:
222  if(fCalculateCorrelations){this->CalculateCorrelations(anEvent);}
223  if(fCalculateDiffCorrelations){this->CalculateDiffCorrelations(anEvent);}
224
225  // f) Calculate e-b-e cumulants: 
226  if(fCalculateEbECumulants){this->CalculateEbECumulants(anEvent);}
227
228  // g) Reset Q-vector components:
229  if(fCalculateQvector||fCalculateDiffQvectors){this->ResetQvector();}
230
231  // h) Cross-check results with nested loops:
232  if(fCrossCheckWithNestedLoops){this->CrossCheckWithNestedLoops(anEvent);}
233  if(fCrossCheckDiffWithNestedLoops){this->CrossCheckDiffWithNestedLoops(anEvent);}
234
235  // i) Dump the points:
236  if(fDumpThePoints){this->DumpThePoints(anEvent);}
237  
238 } // end of AliFlowAnalysisWithMultiparticleCorrelations::Make(AliFlowEventSimple *anEvent)
239
240 //=======================================================================================================================
241
242 void AliFlowAnalysisWithMultiparticleCorrelations::Finish()
243 {
244  // Closing the curtains. 
245
246  // a) Cross-check pointers used in this method;
247  // b) Calculate 'standard candles';
248  // c) Calculate Q-cumulants.
249  
250  // a) Cross-check pointers used in this method:
251  this->CrossCheckPointersUsedInFinish();
252
253  // b) Calculate 'standard candles':
254  if(fCalculateStandardCandles){this->CalculateStandardCandles();}
255
256  // c) Calculate Q-cumulants:
257  if(fCalculateQcumulants){this->CalculateQcumulants();this->CalculateReferenceFlow();}
258
259  // ...
260
261  printf("\n  ... Closing the curtains ... \n\n");
262
263 } // end of AliFlowAnalysisWithMultiparticleCorrelations::Finish()
264
265 //=======================================================================================================================
266
267 void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInFinish()
268 {
269  // Cross-check all pointers used in method Finish().
270
271  // a) Correlations;
272  // b) 'Standard candles';
273  // c) Q-cumulants.
274
275  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInFinish()"; 
276
277  // a) Correlations:
278  if(fCalculateCorrelations)
279  {
280   for(Int_t cs=0;cs<2;cs++) 
281   {
282    if(fCalculateOnlyCos && 1==cs){continue;}
283    else if(fCalculateOnlySin && 0==cs){continue;}
284    for(Int_t c=0;c<fMaxCorrelator;c++) 
285    {
286     if(!fCorrelationsPro[cs][c]){Fatal(sMethodName.Data(),"fCorrelationsPro[%d][%d] is NULL, for one reason or another...",cs,c);}
287    }
288   }
289   if(fCalculateQcumulants && fPropagateErrorQC && !fProductsQCPro)
290    {Fatal(sMethodName.Data(),"fCalculateQcumulants && fPropagateErrorQC && !fProductsQCPro");}
291  } // if(fCalculateCorrelations)
292
293  // b) 'Standard candles':
294  if(fCalculateStandardCandles)
295  {
296   if(!fStandardCandlesHist){Fatal(sMethodName.Data(),"fStandardCandlesHist is NULL, for one reason or another...");}
297   if(fPropagateErrorSC)
298   {
299    if(!fProductsSCPro){Fatal(sMethodName.Data(),"fProductsSCPro is NULL, for one reason or another...");}
300   }
301  } // if(fCalculateStandardCandles)
302
303  // c) Q-cumulants:
304  if(fCalculateQcumulants)
305  {
306   if(!fQcumulantsHist){Fatal(sMethodName.Data(),"fQcumulantsHist is NULL, for one reason or another...");}
307   if(!fReferenceFlowHist){Fatal(sMethodName.Data(),"fReferenceFlowHist is NULL, for one reason or another...");}
308   if(fPropagateErrorQC)
309   {
310    if(!fProductsQCPro){Fatal(sMethodName.Data(),"fProductsQCPro is NULL, for one reason or another...");}
311   }
312  } // if(fCalculateQcumulants)
313
314 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInFinish()
315
316 //=======================================================================================================================
317
318 void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInMake()
319 {
320  // Cross-check all pointers used in method Make().
321
322  // a) Correlations;
323  // b) Event-by-event cumulants;
324  // c) ...
325
326  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInMake()"; 
327
328  // a) Correlations:
329  if(fCalculateCorrelations)
330  {
331   for(Int_t cs=0;cs<2;cs++) 
332   {
333    if(fCalculateOnlyCos && 1==cs){continue;}
334    else if(fCalculateOnlySin && 0==cs){continue;}
335    for(Int_t c=0;c<fMaxCorrelator;c++) 
336    {
337     if(!fCorrelationsPro[cs][c]){Fatal(sMethodName.Data(),"fCorrelationsPro[%d][%d] is NULL, for one reason or another...",cs,c);}
338    }
339   }
340   if(fCalculateQcumulants && fPropagateErrorQC && !fProductsQCPro)
341    {Fatal(sMethodName.Data(),"fCalculateQcumulants && fPropagateErrorQC && !fProductsQCPro");}
342  } // if(fCalculateCorrelations)
343
344  // b) Event-by-event cumulants:
345  if(fCalculateEbECumulants)
346  {
347   for(Int_t cs=0;cs<2;cs++) 
348   {
349    for(Int_t c=0;c<fMaxCorrelator;c++) 
350    {
351     if(!fEbECumulantsPro[cs][c]){Fatal(sMethodName.Data(),"fEbECumulantsPro[%d][%d] is NULL, for one reason or another...",cs,c);}
352    }
353   }
354  } // if(fCalculateEbECumulants)
355
356 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInMake()
357
358 //=======================================================================================================================
359
360 void AliFlowAnalysisWithMultiparticleCorrelations::CalculateStandardCandles()
361 {
362  // Calculate 'standard candles'.
363
364  // 'Standard candle' (SC) is defined in terms of average (all-event!) correlations as follows: 
365  //   SC(-n1,-n2,n2,n1) = <<Cos(-n1,-n2,n2,n1)>> - <<Cos(-n1,n1)>>*<<Cos(-n2,n2)>>, n1 > n2.
366
367  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateStandardCandles()"; 
368
369  Int_t nBins = fStandardCandlesHist->GetNbinsX(); 
370  Int_t nBins2p = fCorrelationsPro[0][1]->GetNbinsX();
371  Int_t nBins4p = fCorrelationsPro[0][3]->GetNbinsX();
372  for(Int_t b=1;b<=nBins;b++)
373  {
374   // Standard candle:
375   Double_t dSCn1n2n2n1 = 0.; // SC(-n1,-n2,n2,n1)
376   Double_t dSCn1n2n2n1Err = 0.; // SC(-n1,-n2,n2,n1) stat. error
377   fPropagateError = kTRUE;
378  
379   // Correlations:
380   Double_t dCosn1n2n2n1 = 0.; // <<Cos(-n1,-n2,n2,n1)>>
381   Double_t dCosn1n2n2n1Err = 0.; // <<Cos(-n1,-n2,n2,n1)>> stat. error
382   Double_t dCosn1n1 = 0.; // <<Cos(-n1,n1)>>
383   Double_t dCosn1n1Err = 0.; // <<Cos(-n1,n1)>> stat. error
384   Double_t dCosn2n2 = 0.; // <<Cos(-n2,n2)>>
385   Double_t dCosn2n2Err = 0.; // <<Cos(-n2,n2)>> stat. error
386   
387   // Labels:
388   TString labeln1n2n2n1 = TString(fStandardCandlesHist->GetXaxis()->GetBinLabel(b)).ReplaceAll("SC","Cos");
389   TString n1 = TString(fStandardCandlesHist->GetXaxis()->GetBinLabel(b))(4);  
390   TString n2 = TString(fStandardCandlesHist->GetXaxis()->GetBinLabel(b))(7);
391   if(n1.EqualTo("-") || n1.EqualTo(",")){Fatal(sMethodName.Data(),"n1.EqualTo...");}
392   if(n2.EqualTo("-") || n2.EqualTo(",")){Fatal(sMethodName.Data(),"n2.EqualTo...");}
393   TString labeln1n1 = Form("Cos(-%s,%s)",n1.Data(),n1.Data());
394   TString labeln2n2 = Form("Cos(-%s,%s)",n2.Data(),n2.Data());
395   //cout<<labeln1n2n2n1.Data()<<endl;
396   //cout<<labeln1n1.Data()<<endl;
397   //cout<<labeln2n2.Data()<<endl;
398   //cout<<endl;  
399
400   // Access <<Cos(-n1,-n2,n2,n1)>>:
401   for(Int_t b4p=1;b4p<=nBins4p;b4p++)
402   {
403    if(labeln1n2n2n1.EqualTo(fCorrelationsPro[0][3]->GetXaxis()->GetBinLabel(b4p)))   
404    {
405     //cout<<labeln1n2n2n1.Data()<<endl;
406     dCosn1n2n2n1 = fCorrelationsPro[0][3]->GetBinContent(b4p);
407     dCosn1n2n2n1Err = fCorrelationsPro[0][3]->GetBinError(b4p);
408     break; 
409    }
410   } // for(Int_t b4p=1;b4p<=nBins4p;b4p++)
411   if(TMath::Abs(dCosn1n2n2n1) < 1.e-44)
412   {
413    cout<<Form("labeln1n2n2n1 = %s",labeln1n2n2n1.Data())<<endl;
414    Warning(sMethodName.Data(),"TMath::Abs(dCosn1n2n2n1) < 1.e-44 !!!!");
415   }
416
417   // Access <<Cos(-n1,n1)>> and <<Cos(-n2,n2)>>:
418   for(Int_t b2p=1;b2p<=nBins2p;b2p++)
419   {
420    if(labeln1n1.EqualTo(fCorrelationsPro[0][1]->GetXaxis()->GetBinLabel(b2p)))   
421    {
422     //cout<<labeln1n1.Data()<<endl;
423     dCosn1n1 = fCorrelationsPro[0][1]->GetBinContent(b2p);
424     dCosn1n1Err = fCorrelationsPro[0][1]->GetBinError(b2p);
425    }  
426    else if(labeln2n2.EqualTo(fCorrelationsPro[0][1]->GetXaxis()->GetBinLabel(b2p)))   
427    {
428     //cout<<labeln2n2.Data()<<endl;
429     dCosn2n2 = fCorrelationsPro[0][1]->GetBinContent(b2p);
430     dCosn2n2Err = fCorrelationsPro[0][1]->GetBinError(b2p);
431    }
432    if(TMath::Abs(dCosn1n1) > 0. && TMath::Abs(dCosn2n2) > 0.){break;} // found 'em both!
433   } // for(Int_t b2p=1;b2p<=nBins2p;b2p++)
434   if(TMath::Abs(dCosn1n1) < 1.e-44)
435   {
436    cout<<Form("labeln1n1 = %s",labeln1n1.Data())<<endl;
437    Warning(sMethodName.Data(),"TMath::Abs(dCosn1n1) < 1.e-44 !!!!");
438   }
439   if(TMath::Abs(dCosn2n2) < 1.e-44)
440   {
441    cout<<Form("labeln2n2 = %s",labeln2n2.Data())<<endl;
442    Warning(sMethodName.Data(),"TMath::Abs(dCosn2n2) < 1.e-44 !!!!");
443   }
444
445   // Calculate standard candles:
446   dSCn1n2n2n1 = dCosn1n2n2n1-dCosn1n1*dCosn2n2;
447
448   // Store the final results:
449   fStandardCandlesHist->SetBinContent(b,dSCn1n2n2n1);
450
451   // Error propagation:
452   if(!fPropagateErrorSC)
453   {
454    fStandardCandlesHist->SetBinError(b,0.);
455    continue;
456   }
457
458   // Access covariances (multiplied by weight dependent prefactor):
459   Double_t wCovCosn1n2n2n1Cosn1n1 = Covariance(labeln1n2n2n1.Data(),labeln1n1.Data(),fProductsSCPro); // weighted Cov(<Cos(-n1,-n2,n2,n1)>,<Cos(-n1,n1)>)
460   Double_t wCovCosn1n2n2n1Cosn2n2 = Covariance(labeln1n2n2n1.Data(),labeln2n2.Data(),fProductsSCPro); // weighted Cov(<Cos(-n1,-n2,n2,n1)>,<Cos(-n2,n2)>)
461   Double_t wCovCosn1n1Cosn2n2 = Covariance(labeln1n1.Data(),labeln2n2.Data(),fProductsSCPro); // weighted Cov(<Cos(-n1,n1)>,<Cos(-n2,n2)>)
462
463   // Explicit error propagation:
464   Double_t dSCn1n2n2n1ErrSquared = pow(dCosn1n1,2.)*pow(dCosn2n2Err,2.) + pow(dCosn2n2,2.)*pow(dCosn1n1Err,2.) 
465                                  + pow(dCosn1n2n2n1Err,2.) + 2.*dCosn1n1*dCosn2n2*wCovCosn1n1Cosn2n2 
466                                  - 2.*dCosn1n1*wCovCosn1n2n2n1Cosn2n2 - 2.*dCosn2n2*wCovCosn1n2n2n1Cosn1n1;
467   if(dSCn1n2n2n1ErrSquared > 0.)
468   {
469    dSCn1n2n2n1Err = pow(dSCn1n2n2n1ErrSquared,0.5);
470   } else
471     {
472      Warning(sMethodName.Data(),"dSCn1n2n2n1ErrSquared > 0. is not satisfied for %s !!!!",labeln1n2n2n1.ReplaceAll("Cos","SC").Data());
473      fPropagateError = kFALSE;
474     }
475
476   // Store the final stat. error:
477   if(fPropagateError)
478   {
479    fStandardCandlesHist->SetBinError(b,dSCn1n2n2n1Err);
480   }
481  } // for(Int_t b=1;b<=nBins;b++)
482
483  fPropagateError = kTRUE;
484
485  return; 
486
487 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateStandardCandles()
488
489 //=======================================================================================================================
490
491 void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForCorrelations()
492 {
493  // Initialize all arrays for correlations.
494
495  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
496  {
497   for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
498   {
499    fCorrelationsPro[cs][c] = NULL;
500   }
501  }
502
503 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForCorrelations()
504
505 //=======================================================================================================================
506
507 void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForEbECumulants()
508 {
509  // Initialize all arrays for event-by-event cumulants.
510
511  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
512  {
513   for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
514   {
515    fEbECumulantsPro[cs][c] = NULL;
516   }
517  }
518
519 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForEbECumulants()
520
521 //=======================================================================================================================
522
523 void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForDiffCorrelations()
524 {
525  // Initialize all arrays for differential correlations.
526
527  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
528  {
529   for(Int_t c=0;c<4;c++) // [1p,2p,3p,4p]
530   {
531    fDiffCorrelationsPro[cs][c] = NULL;
532    fDiffHarmonics[cs][c] = 0;
533   }
534  }
535
536  // Default values:
537  // Cos, 2p:
538  fDiffHarmonics[1][0] = -2;
539  fDiffHarmonics[1][1] = 2;
540  // Cos, 3p:
541  fDiffHarmonics[2][0] = -3;
542  fDiffHarmonics[2][1] = 1;
543  fDiffHarmonics[2][2] = 2;
544  // Cos, 4p:
545  fDiffHarmonics[3][0] = -2;
546  fDiffHarmonics[3][1] = -2;
547  fDiffHarmonics[3][2] = 2;
548  fDiffHarmonics[3][3] = 2;
549
550 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForDiffCorrelations()
551
552 //=======================================================================================================================
553
554 void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForNestedLoops()
555 {
556  // Initialize all arrays for nested loops.  
557
558  fCrossCheckDiffCSCOBN[0] = 0; // cos/sin
559  fCrossCheckDiffCSCOBN[1] = 2; // correlator order
560  fCrossCheckDiffCSCOBN[2] = 4; // bin number
561
562 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForNestedLoops()
563
564 //=======================================================================================================================
565
566 void AliFlowAnalysisWithMultiparticleCorrelations::CalculateCorrelations(AliFlowEventSimple *anEvent)
567 {
568  // Calculate multi-particle correlations from Q-vector components.
569
570  // a) Calculate all booked multi-particle correlations;
571  // b) Calculate products needed for QC error propagation;
572  // c) Calculate products needed for SC error propagation.
573
574  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateCorrelations(AliFlowEventSimple *anEvent)"; 
575  if(!anEvent){Fatal(sMethodName.Data(),"'anEvent'!?!? You again!!!!");}
576
577  // a) Calculate all booked multi-particle correlations:
578  Double_t dMultRP = anEvent->GetNumberOfRPs(); // TBI shall I promote this variable into data member? 
579  for(Int_t cs=0;cs<2;cs++) // cos/sin 
580  {
581   if(fCalculateOnlyCos && 1==cs){continue;}
582   else if(fCalculateOnlySin && 0==cs){continue;}
583   for(Int_t co=0;co<8;co++) // correlator order (TBI hardwired 8) 
584   {
585    if(dMultRP < co+1){break;} // defines min. number of particles in an event for a certain correlator to make sense
586    Int_t nBins = 0;
587    if(fCorrelationsPro[cs][co]){nBins = fCorrelationsPro[cs][co]->GetNbinsX();}
588    else{continue;}
589    for(Int_t b=1;b<=nBins;b++)
590    {
591     TString sBinLabel = fCorrelationsPro[cs][co]->GetXaxis()->GetBinLabel(b);
592     if(sBinLabel.EqualTo("")){break;} 
593     Double_t num = CastStringToCorrelation(sBinLabel.Data(),kTRUE);
594     Double_t den = CastStringToCorrelation(sBinLabel.Data(),kFALSE);
595     Double_t weight = den; // TBI: add support for other options for the weight eventually
596     if(den>0.) 
597     {
598      fCorrelationsPro[cs][co]->Fill(b-.5,num/den,weight);
599     } else{Warning(sMethodName.Data(),"if(den>0.)");}
600    } // for(Int_t b=1;b<=nBins;b++)
601   } // for(Int_t co=0;co<8;co++) // correlator order (TBI hardwired 8) 
602  } // for(Int_t cs=0;cs<=1;cs++) // cos/sin 
603
604  // b) Calculate products needed for QC error propagation:
605  if(fCalculateQcumulants && fPropagateErrorQC){this->CalculateProductsOfCorrelations(anEvent,fProductsQCPro);}
606
607  // c) Calculate products needed for SC error propagation:
608  if(fCalculateStandardCandles && fPropagateErrorSC){this->CalculateProductsOfCorrelations(anEvent,fProductsSCPro);}
609
610 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateCorrelations(AliFlowEventSimple *anEvent)
611
612 //=======================================================================================================================
613
614 void AliFlowAnalysisWithMultiparticleCorrelations::CalculateDiffCorrelations(AliFlowEventSimple *anEvent)
615 {
616  // Calculate differential multi-particle correlations from Q-, p- and q-vector components.
617
618  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateCorrelations(AliFlowEventSimple *anEvent)"; 
619  if(!anEvent){Fatal(sMethodName.Data(),"'anEvent'!?!? You again!!!!");}
620
621  Int_t nBins = 0; // TBI promote this to data member? 
622  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
623  {
624   if(nBins != 0){break;}
625   for(Int_t co=0;co<4;co++) // [1p,2p,3p,4p]
626   {
627    if(fDiffCorrelationsPro[cs][co] && 0==nBins)
628    {
629     nBins = fDiffCorrelationsPro[cs][co]->GetNbinsX(); 
630    }
631   } // for(Int_t co=0;co<4;co++) // [1p,2p,3p,4p]
632  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
633
634  // TBI: The lines below are genuine, most delicious, spaghetti ever... To be reimplemented (one day).
635  if(fCalculateDiffCos)
636  {
637  for(Int_t b=1;b<=nBins;b++)
638  {
639   fDiffBinNo = b-1;
640   // <2'>:  
641   Double_t num2 = TwoDiff(fDiffHarmonics[1][0],fDiffHarmonics[1][1]).Re();
642   Double_t den2 = TwoDiff(0,0).Re();
643   Double_t w2 = den2; // TBI add support for other options for the weight
644   if(den2>0.){fDiffCorrelationsPro[0][1]->Fill(fDiffCorrelationsPro[0][1]->GetBinCenter(b),num2/den2,w2);} 
645   // <3'>:  
646   Double_t num3 = ThreeDiff(fDiffHarmonics[2][0],fDiffHarmonics[2][1],fDiffHarmonics[2][2]).Re();
647   Double_t den3 = ThreeDiff(0,0,0).Re();
648   Double_t w3 = den3; // TBI add support for other options for the weight
649   if(den3>0.){fDiffCorrelationsPro[0][2]->Fill(fDiffCorrelationsPro[0][2]->GetBinCenter(b),num3/den3,w3);} 
650   // <4'>:  
651   Double_t num4 = FourDiff(fDiffHarmonics[3][0],fDiffHarmonics[3][1],fDiffHarmonics[3][2],fDiffHarmonics[3][3]).Re();
652   Double_t den4 = FourDiff(0,0,0,0).Re();
653   Double_t w4 = den4; // TBI add support for other options for the weight
654   if(den4>0.){fDiffCorrelationsPro[0][3]->Fill(fDiffCorrelationsPro[0][3]->GetBinCenter(b),num4/den4,w4);} 
655  } // for(Int_t b=1;b<=nBins;b++)
656  }
657  // TBI: The lines below are genuine, most delicious, spaghetti ever... To be reimplemented (one day).
658  if(fCalculateDiffSin)
659  {
660  for(Int_t b=1;b<=nBins;b++)
661  {
662   fDiffBinNo = b-1;
663   // <2'>:  
664   Double_t num2 = TwoDiff(fDiffHarmonics[1][0],fDiffHarmonics[1][1]).Im();
665   Double_t den2 = TwoDiff(0,0).Re();
666   Double_t w2 = den2; // TBI add support for other options for the weight
667   if(den2>0.){fDiffCorrelationsPro[1][1]->Fill(fDiffCorrelationsPro[1][1]->GetBinCenter(b),num2/den2,w2);} 
668   // <3'>:  
669   Double_t num3 = ThreeDiff(fDiffHarmonics[2][0],fDiffHarmonics[2][1],fDiffHarmonics[2][2]).Im();
670   Double_t den3 = ThreeDiff(0,0,0).Re();
671   Double_t w3 = den3; // TBI add support for other options for the weight
672   if(den3>0.){fDiffCorrelationsPro[1][2]->Fill(fDiffCorrelationsPro[1][2]->GetBinCenter(b),num3/den3,w3);} 
673   // <4'>:  
674   Double_t num4 = FourDiff(fDiffHarmonics[3][0],fDiffHarmonics[3][1],fDiffHarmonics[3][2],fDiffHarmonics[3][3]).Im();
675   Double_t den4 = FourDiff(0,0,0,0).Re();
676   Double_t w4 = den4; // TBI add support for other options for the weight
677   if(den4>0.){fDiffCorrelationsPro[1][3]->Fill(fDiffCorrelationsPro[1][3]->GetBinCenter(b),num4/den4,w4);} 
678  } // for(Int_t b=1;b<=nBins;b++)
679  }
680
681 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateDiffCorrelations(AliFlowEventSimple *anEvent)
682
683 //=======================================================================================================================
684
685 Double_t AliFlowAnalysisWithMultiparticleCorrelations::CastStringToCorrelation(const char *string, Bool_t numerator)
686 {
687  // Cast string of the generic form Cos/Sin(-n_1,-n_2,...,n_{k-1},n_k) in the corresponding correlation value.
688  // If you issue a call to this method with setting numerator = kFALSE, then you are getting back for free
689  // the corresponding denumerator (a.k.a. weight 'number of combinations').
690
691  // TBI:
692  // a) add protection against cases a la:
693  //     string = Cos(-3,-4,5,6,5,6,-3)
694  //     method = Six(-3,-4,5,6,5,-3).Re()
695  // b) cross-check with nested loops this method 
696
697  Double_t dValue = 0.; // return value
698
699  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CastStringToCorrelation(const char *string, Bool_t numerator)"; 
700
701  if(!(TString(string).BeginsWith("Cos") || TString(string).BeginsWith("Sin")))
702  {
703   cout<<Form("And the fatal string is... '%s'. Congratulations!!",string)<<endl; 
704   Fatal(sMethodName.Data(),"!(TString(string).BeginsWith(...");
705  }
706
707  Bool_t bRealPart = kTRUE;
708  if(TString(string).BeginsWith("Sin")){bRealPart = kFALSE;}
709
710  Int_t n[8] = {0,0,0,0,0,0,0,0}; // harmonics, supporting up to 8p correlations
711  UInt_t whichCorr = 0;   
712  for(Int_t t=0;t<=TString(string).Length();t++)
713  {
714   if(TString(string[t]).EqualTo(",") || TString(string[t]).EqualTo(")")) // TBI this is just ugly
715   {
716    n[whichCorr] = string[t-1] - '0';
717    if(TString(string[t-2]).EqualTo("-")){n[whichCorr] = -1*n[whichCorr];}
718    if(!(TString(string[t-2]).EqualTo("-") 
719       || TString(string[t-2]).EqualTo(",")
720       || TString(string[t-2]).EqualTo("("))) // TBI relax this eventually to allow two-digits harmonics
721    { 
722     cout<<Form("And the fatal string is... '%s'. Congratulations!!",string)<<endl; 
723     Fatal(sMethodName.Data(),"!(TString(string[t-2]).EqualTo(...");
724    }
725    whichCorr++;
726    if(whichCorr>=9){Fatal(sMethodName.Data(),"whichCorr>=9");} // not supporting corr. beyond 8p 
727   } // if(TString(string[t]).EqualTo(",") || TString(string[t]).EqualTo(")")) // TBI this is just ugly
728  } // for(UInt_t t=0;t<=TString(string).Length();t++)
729
730  switch(whichCorr)
731  {
732   case 1:
733    if(!numerator){dValue = One(0).Re();}
734    else if(bRealPart){dValue = One(n[0]).Re();} 
735    else{dValue = One(n[0]).Im();}
736   break;
737
738   case 2: 
739    if(!numerator){dValue = Two(0,0).Re();}
740    else if(bRealPart){dValue = Two(n[0],n[1]).Re();}
741    else{dValue = Two(n[0],n[1]).Im();}
742   break;
743
744   case 3: 
745    if(!numerator){dValue = Three(0,0,0).Re();}
746    else if(bRealPart){dValue = Three(n[0],n[1],n[2]).Re();}
747    else{dValue = Three(n[0],n[1],n[2]).Im();}
748   break;
749
750   case 4: 
751    if(!numerator){dValue = Four(0,0,0,0).Re();}
752    else if(bRealPart){dValue = Four(n[0],n[1],n[2],n[3]).Re();}
753    else{dValue = Four(n[0],n[1],n[2],n[3]).Im();}
754   break;
755
756   case 5: 
757    if(!numerator){dValue = Five(0,0,0,0,0).Re();}
758    else if(bRealPart){dValue = Five(n[0],n[1],n[2],n[3],n[4]).Re();}
759    else{dValue = Five(n[0],n[1],n[2],n[3],n[4]).Im();}
760   break;
761
762   case 6: 
763    if(!numerator){dValue = Six(0,0,0,0,0,0).Re();}
764    else if(bRealPart){dValue = Six(n[0],n[1],n[2],n[3],n[4],n[5]).Re();}
765    else{dValue = Six(n[0],n[1],n[2],n[3],n[4],n[5]).Im();}
766   break;
767
768   case 7: 
769    if(!numerator){dValue = Seven(0,0,0,0,0,0,0).Re();}
770    else if(bRealPart){dValue = Seven(n[0],n[1],n[2],n[3],n[4],n[5],n[6]).Re();}
771    else{dValue = Seven(n[0],n[1],n[2],n[3],n[4],n[5],n[6]).Im();}
772    break;
773
774   case 8: 
775    if(!numerator){dValue = Eight(0,0,0,0,0,0,0,0).Re();}
776    else if(bRealPart){dValue = Eight(n[0],n[1],n[2],n[3],n[4],n[5],n[6],n[7]).Re();} 
777    else{dValue = Eight(n[0],n[1],n[2],n[3],n[4],n[5],n[6],n[7]).Im();} 
778   break;
779
780   default:
781    cout<<Form("And the fatal 'whichCorr' value is... %d. Congratulations!!",whichCorr)<<endl; 
782    Fatal(sMethodName.Data(),"switch(whichCorr)"); 
783  } // switch(whichCorr)
784  
785  return dValue;
786
787 } // Double_t AliFlowAnalysisWithMultiparticleCorrelations::CastStringToCorrelation(const char *string, Bool_t numerator)
788
789 //=======================================================================================================================
790
791 void AliFlowAnalysisWithMultiparticleCorrelations::CalculateProductsOfCorrelations(AliFlowEventSimple *anEvent, TProfile2D *profile2D)
792 {
793  // Calculate products of multi-particle correlations (needed for error propagation).
794
795  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateProductsOfCorrelations(AliFlowEventSimple *anEvent, TProfile2D *profile2D)"; 
796  if(!anEvent){Fatal(sMethodName.Data(),"Sorry, 'anEvent' is on holidays.");} 
797  if(!profile2D){Fatal(sMethodName.Data(),"Sorry, 'profile2D' is on holidays.");} 
798
799  Int_t nBins = profile2D->GetXaxis()->GetNbins();
800  for(Int_t bx=2;bx<=nBins;bx++)
801  {
802   for(Int_t by=1;by<bx;by++)
803   {
804    const char *binLabelX = profile2D->GetXaxis()->GetBinLabel(bx);
805    const char *binLabelY = profile2D->GetYaxis()->GetBinLabel(by);
806    Double_t numX = this->CastStringToCorrelation(binLabelX,kTRUE); // numerator
807    Double_t denX = this->CastStringToCorrelation(binLabelX,kFALSE); // denominator
808    Double_t wX = denX; // weight TBI add support for other options
809    Double_t numY = this->CastStringToCorrelation(binLabelY,kTRUE); // numerator
810    Double_t denY = this->CastStringToCorrelation(binLabelY,kFALSE); // denominator
811    Double_t wY = denY; // weight TBI add support for other options
812    if(TMath::Abs(denX) > 0. && TMath::Abs(denY) > 0.)
813    {
814     profile2D->Fill(bx-0.5,by-0.5,(numX/denX)*(numY/denY),wX*wY);
815    } else
816      {
817       cout<<endl; 
818       cout<<"Cannot calculate product for:"<<endl;    
819       cout<<Form("binLabelX = %s",binLabelX)<<endl;
820       cout<<Form("binLabelY = %s",binLabelY)<<endl;
821       cout<<Form("anEvent->GetNumberOfRPs() = %d",anEvent->GetNumberOfRPs())<<endl; 
822       Fatal(sMethodName.Data(),"if(TMath::Abs(denX) > 0. && TMath::Abs(denY) > 0.)");
823      } // else
824   } // for(Int_t by=1;by<bx;by++)
825  } // for(Int_t bx=2;bx<=nBins;bx++)
826
827 } // void CalculateProductsOfCorrelations(AliFlowEventSimple *anEvent, TProfile2D *profile2D)
828
829 //=======================================================================================================================
830
831 void AliFlowAnalysisWithMultiparticleCorrelations::CalculateEbECumulants(AliFlowEventSimple *anEvent)
832 {
833  // Calculate e-b-e cumulants from Q-vector components.
834
835  // TBI this mathod is far (very far, in fact) from being finalized :'(
836
837  // a) Calculate and store e-b-e cumulants.
838
839  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateEbECumulants(AliFlowEventSimple *anEvent)"; 
840  if(!anEvent){Fatal(sMethodName.Data(),"'anEvent'!?!? You again!!!!");}
841
842  // a) Calculate and store e-b-e cumulants:
843  Double_t dMultRP = anEvent->GetNumberOfRPs(); // TBI shall I promote this variable into data member? 
844  Int_t binNo[8]; for(Int_t c=0;c<8;c++){binNo[c]=1;} 
845  // 1-p:
846  for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++) 
847  {
848   if(fSkipZeroHarmonics && 0==n1){continue;}
849   if(fCalculateAll)
850   {
851    TComplex oneN = One(n1); // numerator
852    Double_t oneD = One(0).Re(); // denominator
853    Double_t oneW = oneD; // weight TBI add other possibilities here for the weight
854    if(oneD>0. && dMultRP>=1) 
855    {
856     fEbECumulantsPro[0][0]->Fill(binNo[0]-.5,oneN.Re()/oneD,oneW);
857     fEbECumulantsPro[1][0]->Fill(binNo[0]++-.5,oneN.Im()/oneD,oneW);
858    } else {Warning(sMethodName.Data(),"if(oneD>0. && dMultRP>=1) ");}
859   } 
860   if(1==fDontGoBeyond){continue;}
861   // 2-p:
862   for(Int_t n2=n1;n2<=fMaxHarmonic;n2++) 
863   {
864    if(fSkipZeroHarmonics && 0==n2){continue;}
865    if(fCalculateAll 
866       || (fCalculateIsotropic && 0==n1+n2) 
867       || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2)) 
868       || (fCalculateSameIsotropic && 0==n1+n2 &&  TMath::Abs(n1)==TMath::Abs(n2)))
869    {
870     Double_t cumulants2pCos = Two(n1,n2).Re()/Two(0,0).Re() 
871                             - (One(n1).Re()/One(0).Re())*(One(n2).Re()/One(0).Re())
872                             + (One(n1).Im()/One(0).Re())*(One(n2).Im()/One(0).Re());
873                             
874     Double_t cumulants2pSin = Two(n1,n2).Im()/Two(0,0).Re() 
875                             - (One(n1).Re()/One(0).Re())*(One(n2).Im()/One(0).Re())
876                             - (One(n2).Re()/One(0).Re())*(One(n1).Im()/One(0).Re());
877
878     if(/*twoD>0. &&*/ dMultRP>=2) 
879     {
880      fEbECumulantsPro[0][1]->Fill(binNo[1]-.5,cumulants2pCos,1.);;
881      fEbECumulantsPro[1][1]->Fill(binNo[1]++-.5,cumulants2pSin,1.);;
882     } else {Warning(sMethodName.Data(),"/*twoD>0. &&*/ dMultRP>=2");} 
883    } 
884    if(2==fDontGoBeyond){continue;}
885    
886    /*
887
888    // 3-p:
889    for(Int_t n3=n2;n3<=fMaxHarmonic;n3++) 
890    {
891     if(fSkipZeroHarmonics && 0==n3){continue;}
892     if(fCalculateAll 
893        || (fCalculateIsotropic && 0==n1+n2+n3) 
894        || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3))
895        || (fCalculateSameIsotropic && 0==n1+n2+n3 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)))
896     {
897      TComplex threeN = Three(n1,n2,n3); // numerator
898      Double_t threeD = Three(0,0,0).Re(); // denominator
899      Double_t threeW = threeD; // weight TBI add other possibilities here for the weight
900      if(threeD>0. && dMultRP>=3) 
901      {
902       fEbECumulantsPro[0][2]->Fill(binNo[2]-.5,threeN.Re()/threeD,threeW);
903       fEbECumulantsPro[1][2]->Fill(binNo[2]++-.5,threeN.Im()/threeD,threeW);
904      } else {Warning(sMethodName.Data(),"threeD>0. && dMultRP>=3");} 
905     }
906     if(3==fDontGoBeyond){continue;}
907     // 4-p:
908     for(Int_t n4=n3;n4<=fMaxHarmonic;n4++) 
909     {
910      if(fSkipZeroHarmonics && 0==n4){continue;}
911      if(fCalculateAll 
912         || (fCalculateIsotropic && 0==n1+n2+n3+n4) 
913         || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
914             && TMath::Abs(n1)==TMath::Abs(n4))
915         || (fCalculateSameIsotropic && 0==n1+n2+n3+n4 && TMath::Abs(n1)==TMath::Abs(n2) 
916             && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)))
917      { 
918       TComplex fourN = Four(n1,n2,n3,n4); // numerator
919       Double_t fourD = Four(0,0,0,0).Re(); // denominator
920       Double_t fourW = fourD; // weight TBI add other possibilities here for the weight
921       if(fourD>0. && dMultRP>=4) 
922       {
923        fEbECumulantsPro[0][3]->Fill(binNo[3]-.5,fourN.Re()/fourD,fourW);
924        fEbECumulantsPro[1][3]->Fill(binNo[3]++-.5,fourN.Im()/fourD,fourW);
925       } else {Warning(sMethodName.Data(),"fourD>0. && dMultRP>=4");}
926      }
927      if(4==fDontGoBeyond){continue;}
928      // 5-p:
929      for(Int_t n5=n4;n5<=fMaxHarmonic;n5++) 
930      {
931       if(fSkipZeroHarmonics && 0==n5){continue;}
932       if(fCalculateAll 
933          || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5) 
934          || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
935              && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5))
936          || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5 && TMath::Abs(n1)==TMath::Abs(n2) 
937              && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5)))
938       { 
939        TComplex fiveN = Five(n1,n2,n3,n4,n5); // numerator
940        Double_t fiveD = Five(0,0,0,0,0).Re(); // denominator
941        Double_t fiveW = fiveD; // weight TBI add other possibilities here for the weight
942        if(fiveD>0. && dMultRP>=5) 
943        {
944         fEbECumulantsPro[0][4]->Fill(binNo[4]-.5,fiveN.Re()/fiveD,fiveW);
945         fEbECumulantsPro[1][4]->Fill(binNo[4]++-.5,fiveN.Im()/fiveD,fiveW);
946        } else {Warning(sMethodName.Data(),"fiveD>0. && dMultRP>=5");}
947       } 
948       if(5==fDontGoBeyond){continue;}
949       // 6-p:  
950       for(Int_t n6=n5;n6<=fMaxHarmonic;n6++) 
951       {
952        if(fSkipZeroHarmonics && 0==n6){continue;}
953        if(fCalculateAll 
954           || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6)  
955           || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
956               && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6))
957           || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
958               && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6)))
959        { 
960         TComplex sixN = Six(n1,n2,n3,n4,n5,n6); // numerator
961         Double_t sixD = Six(0,0,0,0,0,0).Re(); // denominator
962         Double_t sixW = sixD; // weight TBI add other possibilities here for the weight
963         if(sixD>0. && dMultRP>=6) 
964         {
965          fEbECumulantsPro[0][5]->Fill(binNo[5]-.5,sixN.Re()/sixD,sixW);
966          fEbECumulantsPro[1][5]->Fill(binNo[5]++-.5,sixN.Im()/sixD,sixW);
967         } else {Warning(sMethodName.Data(),"sixD>0. && dMultRP>=6");}
968        } 
969        if(6==fDontGoBeyond){continue;}
970        // 7-p:
971        for(Int_t n7=n6;n7<=fMaxHarmonic;n7++) 
972        {
973         if(fSkipZeroHarmonics && 0==n7){continue;}
974         if(fCalculateAll 
975            || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7) 
976            || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
977                && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7))
978            || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
979                && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) 
980                && TMath::Abs(n1)==TMath::Abs(n7)))
981         { 
982          TComplex sevenN = Seven(n1,n2,n3,n4,n5,n6,n7); // numerator
983          Double_t sevenD = Seven(0,0,0,0,0,0,0).Re(); // denominator
984          Double_t sevenW = sevenD; // weight TBI add other possibilities here for the weight
985          if(sevenD>0. && dMultRP>=7) 
986          {
987           fEbECumulantsPro[0][6]->Fill(binNo[6]-.5,sevenN.Re()/sevenD,sevenW);
988           fEbECumulantsPro[1][6]->Fill(binNo[6]++-.5,sevenN.Im()/sevenD,sevenW);
989          } else {Warning(sMethodName.Data(),"sevenD>0. && dMultRP>=7");}
990         } 
991         if(7==fDontGoBeyond){continue;}
992         // 8-p:
993         for(Int_t n8=n7;n8<=fMaxHarmonic;n8++) 
994         {
995          if(fSkipZeroHarmonics && 0==n8){continue;}
996          if(fCalculateAll 
997             || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8) 
998             || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
999                 && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7) 
1000                 && TMath::Abs(n1)==TMath::Abs(n8))
1001             || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
1002                 && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) 
1003                 && TMath::Abs(n1)==TMath::Abs(n7) && TMath::Abs(n1)==TMath::Abs(n8)))
1004          { 
1005           TComplex eightN = Eight(n1,n2,n3,n4,n5,n6,n7,n8); // numerator
1006           Double_t eightD = Eight(0,0,0,0,0,0,0,0).Re(); // denominator
1007           Double_t eightW = eightD; // weight TBI add other possibilities here for the weight
1008           if(eightD>0. && dMultRP>=8) 
1009           {
1010            fEbECumulantsPro[0][7]->Fill(binNo[7]-.5,eightN.Re()/eightD,eightW);
1011            fEbECumulantsPro[1][7]->Fill(binNo[7]++-.5,eightN.Im()/eightD,eightW);
1012           }
1013          } 
1014         } // for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
1015        } // for(Int_t n7=n6;n7<=fMaxHarmonic;n7++) 
1016       } // for(Int_t n6=n5;n6<=fMaxHarmonic;n6++) 
1017      } // for(Int_t n5=n4;n5<=fMaxHarmonic;n5++) 
1018     } // for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)   
1019    } // for(Int_t n3=n2;n3<=fMaxHarmonic;n3++) 
1020  
1021   */
1022
1023   } // for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
1024  } // for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++) 
1025
1026 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateEbECumulants(AliFlowEventSimple *anEvent)
1027
1028 //=======================================================================================================================
1029
1030 void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckWithNestedLoops(AliFlowEventSimple *anEvent)
1031 {
1032  // Cross-check results for multi-particle correlations with nested loops.
1033
1034  // TBI add few comments here, there and over there
1035  // TBI this method is rather messy :'(
1036
1037  Int_t h1 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(2);
1038  Int_t h2 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(3);
1039  Int_t h3 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(4);
1040  Int_t h4 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(5);
1041  Int_t h5 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(6);
1042  Int_t h6 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(7);
1043  Int_t h7 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(8);
1044  Int_t h8 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(9);
1045
1046  this->ResetQvector();
1047  this->FillQvector(anEvent);
1048
1049  if(TMath::Abs(One(0).Re())>0.)
1050  {
1051   fNestedLoopsResultsCosPro->Fill(1.5,One(h1).Re()/One(0).Re(),One(0).Re()); 
1052   fNestedLoopsResultsSinPro->Fill(1.5,One(h1).Im()/One(0).Re(),One(0).Re());  
1053  } 
1054  if(TMath::Abs(Two(0,0).Re())>0.)
1055  {
1056   fNestedLoopsResultsCosPro->Fill(3.5,Two(h1,h2).Re()/Two(0,0).Re(),Two(0,0).Re()); 
1057   fNestedLoopsResultsSinPro->Fill(3.5,Two(h1,h2).Im()/Two(0,0).Re(),Two(0,0).Re()); 
1058  }
1059  if(TMath::Abs(Three(0,0,0).Re())>0.)
1060  {
1061   fNestedLoopsResultsCosPro->Fill(5.5,Three(h1,h2,h3).Re()/Three(0,0,0).Re(),Three(0,0,0).Re()); 
1062   fNestedLoopsResultsSinPro->Fill(5.5,Three(h1,h2,h3).Im()/Three(0,0,0).Re(),Three(0,0,0).Re()); 
1063  } 
1064  if(TMath::Abs(Four(0,0,0,0).Re())>0.)
1065  {
1066   fNestedLoopsResultsCosPro->Fill(7.5,Four(h1,h2,h3,h4).Re()/Four(0,0,0,0).Re(),Four(0,0,0,0).Re()); 
1067   fNestedLoopsResultsSinPro->Fill(7.5,Four(h1,h2,h3,h4).Im()/Four(0,0,0,0).Re(),Four(0,0,0,0).Re()); 
1068  } 
1069  if(TMath::Abs(Five(0,0,0,0,0).Re())>0.)
1070  {
1071   fNestedLoopsResultsCosPro->Fill(9.5,Five(h1,h2,h3,h4,h5).Re()/Five(0,0,0,0,0).Re(),Five(0,0,0,0,0).Re()); 
1072   fNestedLoopsResultsSinPro->Fill(9.5,Five(h1,h2,h3,h4,h5).Im()/Five(0,0,0,0,0).Re(),Five(0,0,0,0,0).Re()); 
1073  } 
1074  if(TMath::Abs(Six(0,0,0,0,0,0).Re())>0.)
1075  {
1076   fNestedLoopsResultsCosPro->Fill(11.5,Six(h1,h2,h3,h4,h5,h6).Re()/Six(0,0,0,0,0,0).Re(),Six(0,0,0,0,0,0).Re()); 
1077   fNestedLoopsResultsSinPro->Fill(11.5,Six(h1,h2,h3,h4,h5,h6).Im()/Six(0,0,0,0,0,0).Re(),Six(0,0,0,0,0,0).Re()); 
1078  }
1079  if(TMath::Abs(Seven(0,0,0,0,0,0,0).Re())>0.)
1080  {
1081   fNestedLoopsResultsCosPro->Fill(13.5,Seven(h1,h2,h3,h4,h5,h6,h7).Re()/Seven(0,0,0,0,0,0,0).Re(),Seven(0,0,0,0,0,0,0).Re()); 
1082   fNestedLoopsResultsSinPro->Fill(13.5,Seven(h1,h2,h3,h4,h5,h6,h7).Im()/Seven(0,0,0,0,0,0,0).Re(),Seven(0,0,0,0,0,0,0).Re()); 
1083  }
1084  if(TMath::Abs(Eight(0,0,0,0,0,0,0,0).Re())>0.)
1085  {
1086   fNestedLoopsResultsCosPro->Fill(15.5,Eight(h1,h2,h3,h4,h5,h6,h7,h8).Re()/Eight(0,0,0,0,0,0,0,0).Re(),Eight(0,0,0,0,0,0,0,0).Re()); 
1087   fNestedLoopsResultsSinPro->Fill(15.5,Eight(h1,h2,h3,h4,h5,h6,h7,h8).Im()/Eight(0,0,0,0,0,0,0,0).Re(),Eight(0,0,0,0,0,0,0,0).Re()); 
1088  }
1089
1090  Int_t nPrim = anEvent->NumberOfTracks(); 
1091  Double_t dMultRP = anEvent->GetNumberOfRPs(); // TBI shall I promote this variable into data member? 
1092  AliFlowTrackSimple *aftsTrack = NULL; 
1093  Double_t dPhi1=0.,dPhi2=0.,dPhi3=0.,dPhi4=0.,dPhi5=0.,dPhi6=0.,dPhi7=0.,dPhi8=0.; 
1094  Double_t wPhi1=1.,wPhi2=1.,wPhi3=1.,wPhi4=1.,wPhi5=1.,wPhi6=1.,wPhi7=1.,wPhi8=1.; 
1095
1096  // 1-particle stuff: TBI       
1097  if(dMultRP>=1)
1098  {
1099   for(Int_t i1=0;i1<nPrim;i1++)
1100   {
1101    aftsTrack = anEvent->GetTrack(i1);
1102    if(!(aftsTrack->InRPSelection())){continue;}
1103    dPhi1 = aftsTrack->Phi(); 
1104    if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1105    // Fill:
1106    fNestedLoopsResultsCosPro->Fill(0.5,TMath::Cos(h1*dPhi1),wPhi1); 
1107    fNestedLoopsResultsSinPro->Fill(0.5,TMath::Sin(h1*dPhi1),wPhi1); 
1108   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1109  } // end of if(nPrim>=1) 
1110
1111  // 2-particle correlations:       
1112  if(dMultRP>=2)
1113  {
1114   for(Int_t i1=0;i1<nPrim;i1++)
1115   {
1116    aftsTrack = anEvent->GetTrack(i1);
1117    if(!(aftsTrack->InRPSelection())){continue;}
1118    dPhi1 = aftsTrack->Phi(); 
1119    if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1120    for(Int_t i2=0;i2<nPrim;i2++)
1121    {
1122     if(i2==i1){continue;}
1123     aftsTrack = anEvent->GetTrack(i2);
1124     if(!(aftsTrack->InRPSelection())){continue;}
1125     dPhi2 = aftsTrack->Phi();
1126     if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1127     // Fill:
1128     fNestedLoopsResultsCosPro->Fill(2.5,TMath::Cos(h1*dPhi1+h2*dPhi2),wPhi1*wPhi2); 
1129     fNestedLoopsResultsSinPro->Fill(2.5,TMath::Sin(h1*dPhi1+h2*dPhi2),wPhi1*wPhi2); 
1130    } // end of for(Int_t i2=0;i2<nPrim;i2++)
1131   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1132  } // end of if(nPrim>=2)
1133
1134  // 3-particle correlations:         
1135  if(dMultRP>=3)
1136  {
1137   for(Int_t i1=0;i1<nPrim;i1++)
1138   {
1139    aftsTrack=anEvent->GetTrack(i1);
1140    if(!(aftsTrack->InRPSelection())){continue;}
1141    dPhi1=aftsTrack->Phi();
1142    if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1143    for(Int_t i2=0;i2<nPrim;i2++)
1144    {
1145     if(i2==i1){continue;}
1146     aftsTrack=anEvent->GetTrack(i2);
1147     if(!(aftsTrack->InRPSelection())){continue;}
1148     dPhi2=aftsTrack->Phi();
1149     if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1150     for(Int_t i3=0;i3<nPrim;i3++)
1151     {
1152      if(i3==i1||i3==i2){continue;}
1153      aftsTrack=anEvent->GetTrack(i3);
1154      if(!(aftsTrack->InRPSelection())){continue;}
1155      dPhi3=aftsTrack->Phi();
1156      if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1157      // Fill:
1158      fNestedLoopsResultsCosPro->Fill(4.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3),wPhi1*wPhi2*wPhi3);
1159      fNestedLoopsResultsSinPro->Fill(4.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3),wPhi1*wPhi2*wPhi3);
1160     } // end of for(Int_t i3=0;i3<nPrim;i3++)
1161    } // end of for(Int_t i2=0;i2<nPrim;i2++)
1162   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1163  } // end of if(nPrim>=3)
1164
1165  // 4-particle correlations:
1166  if(dMultRP>=4)
1167  {       
1168   for(Int_t i1=0;i1<nPrim;i1++)
1169   { 
1170    aftsTrack=anEvent->GetTrack(i1);
1171    if(!(aftsTrack->InRPSelection())){continue;}
1172    dPhi1=aftsTrack->Phi();
1173    if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1174    for(Int_t i2=0;i2<nPrim;i2++)
1175    {
1176     if(i2==i1){continue;}
1177     aftsTrack=anEvent->GetTrack(i2);
1178     if(!(aftsTrack->InRPSelection())){continue;}
1179     dPhi2=aftsTrack->Phi();
1180     if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1181     for(Int_t i3=0;i3<nPrim;i3++)
1182     {
1183      if(i3==i1||i3==i2){continue;}
1184      aftsTrack=anEvent->GetTrack(i3);
1185      if(!(aftsTrack->InRPSelection())){continue;}
1186      dPhi3=aftsTrack->Phi();
1187      if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1188      for(Int_t i4=0;i4<nPrim;i4++)
1189      {
1190       if(i4==i1||i4==i2||i4==i3){continue;}
1191       aftsTrack=anEvent->GetTrack(i4);
1192       if(!(aftsTrack->InRPSelection())){continue;}
1193       dPhi4=aftsTrack->Phi();
1194       if(fUseWeights[0][0]){wPhi4 = Weight(dPhi4,"RP","phi");}
1195       // Fill:
1196       fNestedLoopsResultsCosPro->Fill(6.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4),wPhi1*wPhi2*wPhi3*wPhi4);
1197       fNestedLoopsResultsSinPro->Fill(6.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4),wPhi1*wPhi2*wPhi3*wPhi4);
1198      } // end of for(Int_t i4=0;i4<nPrim;i4++) 
1199     } // end of for(Int_t i3=0;i3<nPrim;i3++)
1200    } // end of for(Int_t i2=0;i2<nPrim;i2++)
1201   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1202  } // end of if(nPrim>=)
1203
1204  // 5-particle correlations:      
1205  if(dMultRP>=5)
1206  {
1207   for(Int_t i1=0;i1<nPrim;i1++)
1208   {
1209    aftsTrack=anEvent->GetTrack(i1);
1210    if(!(aftsTrack->InRPSelection())){continue;}  
1211    dPhi1=aftsTrack->Phi();
1212    if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1213    for(Int_t i2=0;i2<nPrim;i2++)
1214    {
1215     if(i2==i1){continue;}
1216     aftsTrack=anEvent->GetTrack(i2);
1217     if(!(aftsTrack->InRPSelection())){continue;}
1218     dPhi2=aftsTrack->Phi();
1219     if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1220     for(Int_t i3=0;i3<nPrim;i3++)
1221     {
1222      if(i3==i1||i3==i2){continue;}
1223      aftsTrack=anEvent->GetTrack(i3);
1224      if(!(aftsTrack->InRPSelection())){continue;}
1225      dPhi3=aftsTrack->Phi();
1226      if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1227      for(Int_t i4=0;i4<nPrim;i4++)
1228      {
1229       if(i4==i1||i4==i2||i4==i3){continue;}
1230       aftsTrack=anEvent->GetTrack(i4);
1231       if(!(aftsTrack->InRPSelection())){continue;}
1232       dPhi4=aftsTrack->Phi();
1233       if(fUseWeights[0][0]){wPhi4 = Weight(dPhi4,"RP","phi");}
1234       for(Int_t i5=0;i5<nPrim;i5++)
1235       {
1236        if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1237        aftsTrack=anEvent->GetTrack(i5);
1238        if(!(aftsTrack->InRPSelection())){continue;}
1239        dPhi5=aftsTrack->Phi();
1240        if(fUseWeights[0][0]){wPhi5 = Weight(dPhi5,"RP","phi");}
1241        // Fill:   
1242        fNestedLoopsResultsCosPro->Fill(8.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5);
1243        fNestedLoopsResultsSinPro->Fill(8.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5);
1244       } // end of for(Int_t i5=0;i5<nPrim;i5++)
1245      } // end of for(Int_t i4=0;i4<nPrim;i4++)  
1246     } // end of for(Int_t i3=0;i3<nPrim;i3++)
1247    } // end of for(Int_t i2=0;i2<nPrim;i2++)
1248   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1249  } // end of if(nPrim>=5)
1250   
1251  // 6-particle correlations:
1252  if(dMultRP>=6)
1253  {
1254   for(Int_t i1=0;i1<nPrim;i1++)
1255   {
1256    aftsTrack=anEvent->GetTrack(i1);
1257    if(!(aftsTrack->InRPSelection())){continue;}
1258    dPhi1=aftsTrack->Phi();
1259    if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1260    for(Int_t i2=0;i2<nPrim;i2++)
1261    {
1262     if(i2==i1){continue;}
1263     aftsTrack=anEvent->GetTrack(i2);
1264     if(!(aftsTrack->InRPSelection())){continue;}
1265     dPhi2=aftsTrack->Phi();
1266     if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1267     for(Int_t i3=0;i3<nPrim;i3++)
1268     {
1269      if(i3==i1||i3==i2){continue;}
1270      aftsTrack=anEvent->GetTrack(i3);
1271      if(!(aftsTrack->InRPSelection())){continue;}
1272      dPhi3=aftsTrack->Phi();
1273      if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1274      for(Int_t i4=0;i4<nPrim;i4++)
1275      {
1276       if(i4==i1||i4==i2||i4==i3){continue;}
1277       aftsTrack=anEvent->GetTrack(i4);
1278       if(!(aftsTrack->InRPSelection())){continue;}
1279       dPhi4=aftsTrack->Phi();
1280       if(fUseWeights[0][0]){wPhi4 = Weight(dPhi4,"RP","phi");}
1281       for(Int_t i5=0;i5<nPrim;i5++)
1282       {
1283        if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1284        aftsTrack=anEvent->GetTrack(i5);
1285        if(!(aftsTrack->InRPSelection())){continue;}
1286        dPhi5=aftsTrack->Phi();
1287        if(fUseWeights[0][0]){wPhi5=Weight(dPhi5,"RP","phi");}
1288        for(Int_t i6=0;i6<nPrim;i6++)
1289        {
1290         if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5){continue;}
1291         aftsTrack=anEvent->GetTrack(i6);
1292         if(!(aftsTrack->InRPSelection())){continue;}
1293         dPhi6=aftsTrack->Phi(); 
1294         if(fUseWeights[0][0]){wPhi6=Weight(dPhi6,"RP","phi");}
1295         // Fill:   
1296         fNestedLoopsResultsCosPro->Fill(10.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6);
1297         fNestedLoopsResultsSinPro->Fill(10.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6);
1298        } // end of for(Int_t i6=0;i6<nPrim;i6++)
1299       } // end of for(Int_t i5=0;i5<nPrim;i5++)
1300      } // end of for(Int_t i4=0;i4<nPrim;i4++)
1301     } // end of for(Int_t i3=0;i3<nPrim;i3++)
1302    } // end of for(Int_t i2=0;i2<nPrim;i2++)
1303   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1304  } // end of if(nPrim>=6)
1305   
1306  // 7-particle correlations:
1307  if(dMultRP>=7)
1308  {
1309   for(Int_t i1=0;i1<nPrim;i1++)
1310   { 
1311    aftsTrack=anEvent->GetTrack(i1);
1312    if(!(aftsTrack->InRPSelection())){continue;}
1313    dPhi1=aftsTrack->Phi();
1314    if(fUseWeights[0][0]){wPhi1=Weight(dPhi1,"RP","phi");}
1315    for(Int_t i2=0;i2<nPrim;i2++)
1316    {
1317     if(i2==i1){continue;}
1318     aftsTrack=anEvent->GetTrack(i2);
1319     if(!(aftsTrack->InRPSelection())){continue;}
1320     dPhi2=aftsTrack->Phi();
1321     if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1322     for(Int_t i3=0;i3<nPrim;i3++)
1323     {
1324      if(i3==i1||i3==i2){continue;}
1325      aftsTrack=anEvent->GetTrack(i3);
1326      if(!(aftsTrack->InRPSelection())){continue;}
1327      dPhi3=aftsTrack->Phi();
1328      if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1329      for(Int_t i4=0;i4<nPrim;i4++)
1330      {
1331       if(i4==i1||i4==i2||i4==i3){continue;}
1332       aftsTrack=anEvent->GetTrack(i4);
1333       if(!(aftsTrack->InRPSelection())){continue;}
1334       dPhi4=aftsTrack->Phi();
1335       if(fUseWeights[0][0]){wPhi4=Weight(dPhi4,"RP","phi");}
1336       for(Int_t i5=0;i5<nPrim;i5++)
1337       {
1338        if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1339        aftsTrack=anEvent->GetTrack(i5);
1340        if(!(aftsTrack->InRPSelection())){continue;}
1341        dPhi5=aftsTrack->Phi();
1342        if(fUseWeights[0][0]){wPhi5=Weight(dPhi5,"RP","phi");}
1343        for(Int_t i6=0;i6<nPrim;i6++)
1344        {
1345         if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5){continue;}
1346         aftsTrack=anEvent->GetTrack(i6);
1347         if(!(aftsTrack->InRPSelection())){continue;}
1348         dPhi6=aftsTrack->Phi(); 
1349         if(fUseWeights[0][0]){wPhi6=Weight(dPhi6,"RP","phi");}
1350         for(Int_t i7=0;i7<nPrim;i7++)
1351         {
1352          if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6){continue;}
1353          aftsTrack=anEvent->GetTrack(i7);
1354          if(!(aftsTrack->InRPSelection())){continue;}
1355          dPhi7=aftsTrack->Phi(); 
1356          if(fUseWeights[0][0]){wPhi7=Weight(dPhi7,"RP","phi");}
1357          // Fill:   
1358          fNestedLoopsResultsCosPro->Fill(12.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6+h7*dPhi7),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6*wPhi7);
1359          fNestedLoopsResultsSinPro->Fill(12.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6+h7*dPhi7),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6*wPhi7);
1360         } // end of for(Int_t i7=0;i7<nPrim;i7++)
1361        } // end of for(Int_t i6=0;i6<nPrim;i6++) 
1362       } // end of for(Int_t i5=0;i5<nPrim;i5++)
1363      } // end of for(Int_t i4=0;i4<nPrim;i4++)  
1364     } // end of for(Int_t i3=0;i3<nPrim;i3++)
1365    } // end of for(Int_t i2=0;i2<nPrim;i2++)
1366   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1367  } // end of if(nPrim>=7)
1368  
1369  // 8-particle correlations:
1370  if(dMultRP>=8)
1371  {
1372   for(Int_t i1=0;i1<nPrim;i1++)
1373   {
1374    aftsTrack=anEvent->GetTrack(i1);
1375    if(!(aftsTrack->InRPSelection())){continue;}
1376    dPhi1=aftsTrack->Phi();
1377    if(fUseWeights[0][0]){wPhi1=Weight(dPhi1,"RP","phi");}
1378    for(Int_t i2=0;i2<nPrim;i2++)
1379    {
1380     if(i2==i1){continue;}
1381     aftsTrack=anEvent->GetTrack(i2);
1382     if(!(aftsTrack->InRPSelection())){continue;}
1383     dPhi2=aftsTrack->Phi();
1384     if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1385     for(Int_t i3=0;i3<nPrim;i3++)
1386     {
1387      if(i3==i1||i3==i2){continue;}
1388      aftsTrack=anEvent->GetTrack(i3);
1389      if(!(aftsTrack->InRPSelection())){continue;}
1390      dPhi3=aftsTrack->Phi();
1391      if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1392      for(Int_t i4=0;i4<nPrim;i4++)
1393      {
1394       if(i4==i1||i4==i2||i4==i3){continue;}
1395       aftsTrack=anEvent->GetTrack(i4);
1396       if(!(aftsTrack->InRPSelection())){continue;}
1397       dPhi4=aftsTrack->Phi();
1398       if(fUseWeights[0][0]){wPhi4=Weight(dPhi4,"RP","phi");}
1399       for(Int_t i5=0;i5<nPrim;i5++)
1400       {
1401        if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1402        aftsTrack=anEvent->GetTrack(i5);
1403        if(!(aftsTrack->InRPSelection())){continue;}
1404        dPhi5=aftsTrack->Phi();
1405        if(fUseWeights[0][0]){wPhi5=Weight(dPhi5,"RP","phi");}
1406        for(Int_t i6=0;i6<nPrim;i6++)
1407        {
1408         if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5){continue;}
1409         aftsTrack=anEvent->GetTrack(i6);
1410         if(!(aftsTrack->InRPSelection())){continue;}
1411         dPhi6=aftsTrack->Phi();
1412         if(fUseWeights[0][0]){wPhi6=Weight(dPhi6,"RP","phi");}
1413         for(Int_t i7=0;i7<nPrim;i7++)
1414         {
1415          if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6){continue;}
1416          aftsTrack=anEvent->GetTrack(i7);
1417          if(!(aftsTrack->InRPSelection())){continue;}
1418          dPhi7=aftsTrack->Phi();
1419          if(fUseWeights[0][0]){wPhi7=Weight(dPhi7,"RP","phi");}
1420          for(Int_t i8=0;i8<nPrim;i8++)
1421          {
1422           if(i8==i1||i8==i2||i8==i3||i8==i4||i8==i5||i8==i6||i8==i7){continue;}
1423           aftsTrack=anEvent->GetTrack(i8);
1424           if(!(aftsTrack->InRPSelection())){continue;}
1425           dPhi8=aftsTrack->Phi();
1426           if(fUseWeights[0][0]){wPhi8=Weight(dPhi8,"RP","phi");}
1427           // Fill:   
1428           fNestedLoopsResultsCosPro->Fill(14.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6+h7*dPhi7+h8*dPhi8),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6*wPhi7*wPhi8);
1429           fNestedLoopsResultsSinPro->Fill(14.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6+h7*dPhi7+h8*dPhi8),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6*wPhi7*wPhi8);
1430          } // end of for(Int_t i8=0;i8<nPrim;i8++)
1431         } // end of for(Int_t i7=0;i7<nPrim;i7++) 
1432        } // end of for(Int_t i6=0;i6<nPrim;i6++) 
1433       } // end of for(Int_t i5=0;i5<nPrim;i5++)
1434      } // end of for(Int_t i4=0;i4<nPrim;i4++)  
1435     } // end of for(Int_t i3=0;i3<nPrim;i3++)
1436    } // end of for(Int_t i2=0;i2<nPrim;i2++)
1437   } // end of for(Int_t i1=0;i1<nPrim;i1++)
1438  } // end of if(nPrim>=8)
1439  
1440  // *) Printout: TBI move somewhere else
1441  printf("\n cosine:");
1442  printf("\n  1-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(2));
1443  printf("\n  1-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(1));
1444  printf("\n  2-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(4));
1445  printf("\n  2-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(3));
1446  printf("\n  3-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(6));
1447  printf("\n  3-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(5));
1448  printf("\n  4-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(8));
1449  printf("\n  4-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(7));
1450  printf("\n  5-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(10));
1451  printf("\n  5-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(9));
1452  printf("\n  6-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(12));
1453  printf("\n  6-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(11));
1454  printf("\n  7-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(14));
1455  printf("\n  7-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(13));
1456  printf("\n  8-p => Q-vector:     %.12f",fNestedLoopsResultsCosPro->GetBinContent(16));
1457  printf("\n  8-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(15));
1458
1459  printf("\n\n sinus:");
1460  printf("\n  1-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(2));
1461  printf("\n  1-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(1));
1462  printf("\n  2-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(4));
1463  printf("\n  2-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(3));
1464  printf("\n  3-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(6));
1465  printf("\n  3-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(5));
1466  printf("\n  4-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(8));
1467  printf("\n  4-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(7));
1468  printf("\n  5-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(10));
1469  printf("\n  5-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(9));
1470  printf("\n  6-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(12));
1471  printf("\n  6-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(11));
1472  printf("\n  7-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(14));
1473  printf("\n  7-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(13));
1474  printf("\n  8-p => Q-vector:     %.12f",fNestedLoopsResultsSinPro->GetBinContent(16));
1475  printf("\n  8-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(15));
1476
1477  printf("\n\n"); 
1478
1479 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckWithNestedLoops(AliFlowEventSimple *anEvent)
1480
1481 //=======================================================================================================================
1482
1483 void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckDiffWithNestedLoops(AliFlowEventSimple *anEvent)
1484 {
1485  // Cross-check results for differential multi-particle correlations with nested loops.
1486
1487  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckDiffWithNestedLoops(AliFlowEventSimple *anEvent)";
1488
1489  Int_t nPrim = anEvent->NumberOfTracks(); 
1490  AliFlowTrackSimple *aftsTrack = NULL; 
1491  Double_t dPsi1=0.,dPhi2=0.,dPhi3=0.,dPhi4=0.; 
1492  Double_t wPsi1=1.,wPhi2=1.,wPhi3=1.,wPhi4=1.; 
1493
1494  Int_t cs = fCrossCheckDiffCSCOBN[0]; // cos/sin
1495
1496  // TBI reimplement lines below in a more civilised manner:
1497  Bool_t bCrossCheck2p = kFALSE;
1498  Bool_t bCrossCheck3p = kFALSE;
1499  Bool_t bCrossCheck4p = kFALSE;
1500
1501  if(fCrossCheckDiffCSCOBN[1] == 2){bCrossCheck2p = kTRUE;}
1502  else if(fCrossCheckDiffCSCOBN[1] == 3){bCrossCheck3p = kTRUE;}
1503  else if(fCrossCheckDiffCSCOBN[1] == 4){bCrossCheck4p = kTRUE;}
1504
1505  if(Int_t(bCrossCheck2p + bCrossCheck3p + bCrossCheck4p) > 1)
1506  {
1507   Fatal(sMethodName.Data(),"Int_t(bCrossCheck2p + bCrossCheck3p + bCrossCheck4p) > 1");
1508  }
1509  if(!(bCrossCheck2p || bCrossCheck3p || bCrossCheck4p))
1510  {
1511   Fatal(sMethodName.Data(),"!(bCrossCheck2p || bCrossCheck3p || bCrossCheck4p)");
1512  }
1513  Int_t nDiffBinNo = fCrossCheckDiffCSCOBN[2];
1514  Double_t dPt = 0., dEta = 0.;
1515
1516  // <2'>:
1517  for(Int_t i1=0;i1<nPrim;i1++) // Loop over particles in a differential bin 
1518  {
1519   aftsTrack=anEvent->GetTrack(i1);
1520   if(!(aftsTrack->InPOISelection())){continue;}
1521   dPsi1=aftsTrack->Phi();
1522   if(fCalculateDiffCorrelationsVsPt)
1523   {
1524    dPt=aftsTrack->Pt();
1525    if(fDiffCorrelationsPro[0][1]->FindBin(dPt) != nDiffBinNo){continue;} // TBI spaghetti again 
1526   } else 
1527     {
1528      dEta=aftsTrack->Eta();
1529      if(fDiffCorrelationsPro[0][1]->FindBin(dEta) != nDiffBinNo){continue;} // TBI spaghetti again 
1530     }
1531   if(fUseWeights[1][0]){wPsi1=Weight(dPsi1,"POI","phi");}
1532   for(Int_t i2=0;i2<nPrim;i2++) // Loop over particles in an event
1533   {
1534    if(i2==i1){continue;} // get rid of autocorrelations
1535    aftsTrack=anEvent->GetTrack(i2);
1536    if(!(aftsTrack->InRPSelection())){continue;}
1537    dPhi2=aftsTrack->Phi();
1538    if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1539    // Fill profiles:
1540    if(bCrossCheck2p)
1541    {
1542     if(fCrossCheckDiffCSCOBN[0] == 0)
1543     {
1544      fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Cos(fDiffHarmonics[1][0]*dPsi1+fDiffHarmonics[1][1]*dPhi2),wPsi1*wPhi2);
1545     } else {fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Sin(fDiffHarmonics[1][0]*dPsi1+fDiffHarmonics[1][1]*dPhi2),wPsi1*wPhi2);}
1546    } // if(bCrossCheck2p) 
1547   } // for(Int_t i2=0;i2<nPrim;i2++)
1548  } // for(Int_t i1=0;i1<nPrim;i1++)
1549
1550  // <3'>:
1551  for(Int_t i1=0;i1<nPrim;i1++) // Loop over particles in a differential bin
1552  {
1553   aftsTrack=anEvent->GetTrack(i1);
1554   if(!(aftsTrack->InPOISelection())){continue;}
1555   dPsi1=aftsTrack->Phi();
1556   if(fCalculateDiffCorrelationsVsPt)
1557   {
1558    dPt=aftsTrack->Pt();
1559    if(fDiffCorrelationsPro[0][1]->FindBin(dPt) != nDiffBinNo){continue;} // TBI spaghetti again 
1560   } else 
1561     {
1562      dEta=aftsTrack->Eta();
1563      if(fDiffCorrelationsPro[0][1]->FindBin(dEta) != nDiffBinNo){continue;} // TBI spaghetti again 
1564     }
1565   if(fUseWeights[1][0]){wPsi1=Weight(dPsi1,"POI","phi");}
1566   for(Int_t i2=0;i2<nPrim;i2++) // Loop over particles in an event
1567   {
1568    if(i2==i1){continue;} // get rid of autocorrelations
1569    aftsTrack=anEvent->GetTrack(i2);
1570    if(!(aftsTrack->InRPSelection())){continue;}
1571    dPhi2=aftsTrack->Phi();
1572    if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1573    for(Int_t i3=0;i3<nPrim;i3++)
1574    {
1575     if(i3==i1||i3==i2){continue;} // get rid of autocorrelations
1576     aftsTrack=anEvent->GetTrack(i3);
1577     if(!(aftsTrack->InRPSelection())){continue;}
1578     dPhi3=aftsTrack->Phi();
1579     if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1580     // Fill the profiles:
1581     if(bCrossCheck3p)
1582     {
1583      if(fCrossCheckDiffCSCOBN[0] == 0)
1584      {
1585       fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Cos(fDiffHarmonics[2][0]*dPsi1+fDiffHarmonics[2][1]*dPhi2+fDiffHarmonics[2][2]*dPhi3),wPsi1*wPhi2*wPhi3);  
1586      } else {fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Sin(fDiffHarmonics[2][0]*dPsi1+fDiffHarmonics[2][1]*dPhi2+fDiffHarmonics[2][2]*dPhi3),wPsi1*wPhi2*wPhi3);}
1587     } // if(bCrossCheck3p)
1588    } // end of for(Int_t i3=0;i3<nPrim;i3++)  
1589   } // for(Int_t i2=0;i2<nPrim;i2++)
1590  } // for(Int_t i1=0;i1<nPrim;i1++)
1591
1592  // <4'>:
1593  for(Int_t i1=0;i1<nPrim;i1++) // Loop over particles in a differential bin
1594  {
1595   aftsTrack=anEvent->GetTrack(i1);
1596   if(!(aftsTrack->InPOISelection())){continue;}
1597   dPsi1=aftsTrack->Phi();
1598   if(fCalculateDiffCorrelationsVsPt)
1599   {
1600    dPt=aftsTrack->Pt();
1601    if(fDiffCorrelationsPro[0][1]->FindBin(dPt) != nDiffBinNo){continue;} // TBI spaghetti again 
1602   } else 
1603     {
1604      dEta=aftsTrack->Eta();
1605      if(fDiffCorrelationsPro[0][1]->FindBin(dEta) != nDiffBinNo){continue;} // TBI spaghetti again 
1606     }
1607   if(fUseWeights[1][0]){wPsi1=Weight(dPsi1,"POI","phi");}
1608   for(Int_t i2=0;i2<nPrim;i2++) // Loop over particles in an event
1609   {
1610    if(i2==i1){continue;} // get rid of autocorrelations
1611    aftsTrack=anEvent->GetTrack(i2);
1612    if(!(aftsTrack->InRPSelection())){continue;}
1613    dPhi2=aftsTrack->Phi();
1614    if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1615    for(Int_t i3=0;i3<nPrim;i3++)
1616    {
1617     if(i3==i1||i3==i2){continue;} // get rid of autocorrelations
1618     aftsTrack=anEvent->GetTrack(i3);
1619     if(!(aftsTrack->InRPSelection())){continue;}
1620     dPhi3=aftsTrack->Phi();
1621     if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1622     for(Int_t i4=0;i4<nPrim;i4++)
1623     {
1624      if(i4==i1||i4==i2||i4==i3){continue;} // get rid of autocorrelations
1625      aftsTrack=anEvent->GetTrack(i4);
1626      if(!(aftsTrack->InRPSelection())){continue;}
1627      dPhi4=aftsTrack->Phi();
1628      if(fUseWeights[0][0]){wPhi4=Weight(dPhi4,"RP","phi");}
1629      // Fill the profiles:
1630      if(bCrossCheck4p)
1631      {
1632       if(fCrossCheckDiffCSCOBN[0] == 0)
1633       {
1634        fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Cos(fDiffHarmonics[3][0]*dPsi1+fDiffHarmonics[3][1]*dPhi2+fDiffHarmonics[3][2]*dPhi3+fDiffHarmonics[3][3]*dPhi4),wPsi1*wPhi2*wPhi3*wPhi4);
1635       } else {fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Sin(fDiffHarmonics[3][0]*dPsi1+fDiffHarmonics[3][1]*dPhi2+fDiffHarmonics[3][2]*dPhi3+fDiffHarmonics[3][3]*dPhi4),wPsi1*wPhi2*wPhi3*wPhi4);} 
1636      } // if(bCrossCheck4p)
1637     } // end of for(Int_t i4=0;i4<nPrim;i4++) 
1638    } // end of for(Int_t i3=0;i3<nPrim;i3++)  
1639   } // for(Int_t i2=0;i2<nPrim;i2++)
1640  } // for(Int_t i1=0;i1<nPrim;i1++)
1641
1642  // Printout:
1643  // 2-p:
1644  if(bCrossCheck2p)
1645  {
1646   printf("\n  2-p => Q-vector:     %.12f",fDiffCorrelationsPro[cs][1]->GetBinContent(nDiffBinNo));
1647   printf("\n  2-p => Nested loops: %.12f\n",fNestedLoopsDiffResultsPro->GetBinContent(1));
1648  }
1649  // 3-p:
1650  if(bCrossCheck3p)
1651  {
1652   printf("\n  3-p => Q-vector:     %.12f",fDiffCorrelationsPro[cs][2]->GetBinContent(nDiffBinNo));
1653   printf("\n  3-p => Nested loops: %.12f\n",fNestedLoopsDiffResultsPro->GetBinContent(1));
1654  } 
1655  // 4-p:
1656  if(bCrossCheck4p)
1657  {
1658   printf("\n  4-p => Q-vector:     %.12f",fDiffCorrelationsPro[cs][3]->GetBinContent(nDiffBinNo));
1659   printf("\n  4-p => Nested loops: %.12f\n",fNestedLoopsDiffResultsPro->GetBinContent(1));
1660  }
1661
1662 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckDiffWithNestedLoops(AliFlowEventSimple *anEvent)
1663
1664 //=======================================================================================================================
1665
1666 void AliFlowAnalysisWithMultiparticleCorrelations::FillQvector(AliFlowEventSimple *anEvent)
1667 {
1668  // Fill Q-vector components.
1669
1670  Int_t nTracks = anEvent->NumberOfTracks(); // TBI shall I promote this to data member?
1671  Double_t dPhi = 0., wPhi = 1.; // azimuthal angle and corresponding phi weight
1672  Double_t dPt = 0., wPt = 1.; // transverse momentum and corresponding pT weight
1673  Double_t dEta = 0., wEta = 1.; // pseudorapidity and corresponding eta weight
1674  Double_t wToPowerP = 1.; // weight raised to power p
1675  for(Int_t t=0;t<nTracks;t++) // loop over all tracks
1676  {
1677   AliFlowTrackSimple *pTrack = anEvent->GetTrack(t);
1678   if(!pTrack){printf("\n AAAARGH: pTrack is NULL in MPC::FillQvector(...) !!!!"); continue;}
1679   if(!(pTrack->InRPSelection() || pTrack->InPOISelection())){printf("\n AAAARGH: pTrack is neither RP nor POI !!!!"); continue;}
1680   if(pTrack->InRPSelection()) // fill Q-vector components only with reference particles
1681   {
1682    wPhi = 1.; wPt = 1.; wEta = 1.; wToPowerP = 1.; // TBI this shall go somewhere else, for performance sake
1683
1684    // Access kinematic variables for RP and corresponding weights:
1685    dPhi = pTrack->Phi(); // azimuthal angle
1686    if(fUseWeights[0][0]){wPhi = Weight(dPhi,"RP","phi");} // corresponding phi weight
1687    //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
1688    //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
1689    dPt = pTrack->Pt();
1690    if(fUseWeights[0][1]){wPt = Weight(dPt,"RP","pt");} // corresponding pT weight
1691    dEta = pTrack->Eta();
1692    if(fUseWeights[0][2]){wEta = Weight(dEta,"RP","eta");} // corresponding eta weight
1693    // Calculate Q-vector components:
1694    for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
1695    {
1696     for(Int_t wp=0;wp<fMaxCorrelator+1;wp++) // weight power
1697     {
1698      if(fUseWeights[0][0]||fUseWeights[0][1]||fUseWeights[0][2]){wToPowerP = pow(wPhi*wPt*wEta,wp);} 
1699      fQvector[h][wp] += TComplex(wToPowerP*TMath::Cos(h*dPhi),wToPowerP*TMath::Sin(h*dPhi));
1700     } // for(Int_t wp=0;wp<fMaxCorrelator+1;wp++)
1701    } // for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
1702   } // if(pTrack->InRPSelection()) // fill Q-vector components only with reference particles
1703
1704   // Differential Q-vectors (a.k.a. p-vector and q-vector):
1705   if(!fCalculateDiffQvectors){continue;}
1706   if(pTrack->InPOISelection()) 
1707   {
1708    wPhi = 1.; wPt = 1.; wEta = 1.; wToPowerP = 1.; // TBI this shall go somewhere else, for performance sake
1709
1710    // Access kinematic variables for POI and corresponding weights:
1711    dPhi = pTrack->Phi(); // azimuthal angle
1712    if(fUseWeights[1][0]){wPhi = Weight(dPhi,"POI","phi");} // corresponding phi weight
1713    //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
1714    //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
1715    dPt = pTrack->Pt();
1716    if(fUseWeights[1][1]){wPt = Weight(dPt,"POI","pt");} // corresponding pT weight
1717    dEta = pTrack->Eta();
1718    if(fUseWeights[1][2]){wEta = Weight(dEta,"POI","eta");} // corresponding eta weight
1719    // Determine bin:
1720    Int_t binNo = -44;
1721    if(fCalculateDiffCorrelationsVsPt)
1722    { 
1723     binNo = fDiffCorrelationsPro[0][0]->FindBin(dPt); // TBI: hardwired [0][0]
1724    } else
1725      {
1726       binNo = fDiffCorrelationsPro[0][0]->FindBin(dEta); // TBI: hardwired [0][0]
1727      }
1728    // Calculate p-vector components:
1729    for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
1730    {
1731     for(Int_t wp=0;wp<fMaxCorrelator+1;wp++) // weight power
1732     {
1733      if(fUseWeights[1][0]||fUseWeights[1][1]||fUseWeights[1][2]){wToPowerP = pow(wPhi*wPt*wEta,wp);} 
1734      fpvector[binNo-1][h][wp] += TComplex(wToPowerP*TMath::Cos(h*dPhi),wToPowerP*TMath::Sin(h*dPhi));
1735
1736      if(pTrack->InRPSelection()) 
1737      {
1738       // Fill q-vector components:
1739       wPhi = 1.; wPt = 1.; wEta = 1.; wToPowerP = 1.; // TBI this shall go somewhere else, for performance sake
1740
1741       if(fUseWeights[0][0]){wPhi = Weight(dPhi,"RP","phi");} // corresponding phi weight
1742       //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
1743       //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
1744       if(fUseWeights[0][1]){wPt = Weight(dPt,"RP","pt");} // corresponding pT weight
1745       if(fUseWeights[0][2]){wEta = Weight(dEta,"RP","eta");} // corresponding eta weight
1746       if(fUseWeights[1][0]){wPhi = Weight(dPhi,"POI","phi");} // corresponding phi weight
1747       //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
1748       //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
1749       if(fUseWeights[1][1]){wPt = Weight(dPt,"POI","pt");} // corresponding pT weight
1750       if(fUseWeights[1][2]){wEta = Weight(dEta,"POI","eta");} // corresponding eta weight
1751       if(fUseWeights[0][0]||fUseWeights[0][1]||fUseWeights[0][2]||fUseWeights[1][0]||fUseWeights[1][1]||fUseWeights[1][2]){wToPowerP = pow(wPhi*wPt*wEta,wp);} 
1752       fqvector[binNo-1][h][wp] += TComplex(wToPowerP*TMath::Cos(h*dPhi),wToPowerP*TMath::Sin(h*dPhi));
1753      } // if(pTrack->InRPSelection()) 
1754
1755     } // for(Int_t wp=0;wp<fMaxCorrelator+1;wp++)
1756    } // for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
1757   } // if(pTrack->InPOISelection()) 
1758
1759  } // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
1760
1761 } // void AliFlowAnalysisWithMultiparticleCorrelations::FillQvector(AliFlowEventSimple *anEvent)
1762
1763 //=======================================================================================================================
1764
1765 void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckSettings()
1766 {
1767  // Cross-check all initial settings in this method. 
1768  
1769  // a) Few cross-checks for control histograms;
1770  // b) Few cross-checks for flags for correlations;
1771  // c) 'Standard candles';
1772  // d) Q-cumulants;
1773  // e) Weights;
1774  // f) Differential correlations;
1775  // g) Nested loops;
1776  // h) Dump the points.
1777
1778  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckSettings()";
1779
1780  // a) Few cross-checks for control histograms: TBI the lines below are not really what they are supposed to be...
1781  /*
1782  if(fFillKinematicsHist && !fFillControlHistograms){Fatal(sMethodName.Data(),"fFillKinematicsHist && !fFillControlHistograms");}
1783  if(fFillMultDistributionsHist && !fFillControlHistograms){Fatal(sMethodName.Data(),"fFillMultDistributionsHist && !fFillControlHistograms");}
1784  if(fFillMultCorrelationsHist && !fFillControlHistograms){Fatal(sMethodName.Data(),"fFillMultCorrelationsHist && !fFillControlHistograms");}
1785  */ 
1786
1787  // b) Few cross-checks for flags for correlations: // TBI the lines bellow can be civilized
1788  Int_t iSum = (Int_t)fCalculateIsotropic + (Int_t)fCalculateSame + (Int_t)fCalculateSameIsotropic;
1789  if(iSum>1){Fatal(sMethodName.Data(),"iSum is doing crazy things...");}
1790  if(fCalculateOnlyCos && fCalculateOnlySin){Fatal(sMethodName.Data(),"fCalculateOnlyCos && fCalculateOnlySin");}
1791
1792  // c) 'Standard candles':
1793  if(fCalculateStandardCandles && !fCalculateCorrelations)
1794  {
1795   Fatal(sMethodName.Data(),"fCalculateStandardCandles && !fCalculateCorrelations");
1796  }
1797  if(fCalculateStandardCandles && fCalculateCorrelations && fCalculateSameIsotropic)
1798  {
1799   Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateCorrelations && fCalculateSameIsotropic");
1800  }
1801  if(fCalculateStandardCandles && fCalculateOnlyForHarmonicQC)
1802  {
1803   Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateOnlyForHarmonicQC");
1804  }
1805  if(fCalculateStandardCandles && fCalculateOnlyForSC && (4!=fDontGoBeyond))
1806  {
1807   Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateOnlyForSC && (4!=fDontGoBeyond)");
1808  }
1809  if(fCalculateStandardCandles && !fPropagateErrorSC)
1810  {
1811   Warning(sMethodName.Data(),"fCalculateStandardCandles && !fPropagateErrorSC");
1812  }
1813  if(fCalculateStandardCandles && fCalculateOnlySin)
1814  {
1815   Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateOnlySin");
1816  }
1817  if(fCalculateStandardCandles && fDontGoBeyond < 3)
1818  {
1819   Fatal(sMethodName.Data(),"fCalculateStandardCandles && fDontGoBeyond < 3");
1820  }
1821
1822  // d) Q-cumulants:
1823  if(fCalculateQcumulants && !fCalculateCorrelations)
1824  {
1825   Fatal(sMethodName.Data(),"fCalculateQcumulants && !fCalculateCorrelations");
1826  }
1827  if(fCalculateQcumulants && !(fHarmonicQC > 0))
1828  {
1829   Fatal(sMethodName.Data(),"fCalculateQcumulants && !(fHarmonicQC > 0)");
1830  }
1831  if(fCalculateQcumulants && fCalculateOnlyForSC)
1832  {
1833   Fatal(sMethodName.Data(),"fCalculateQcumulants && fCalculateOnlyForSC");
1834  }
1835  if(fCalculateQcumulants && !fPropagateErrorQC)
1836  {
1837   Warning(sMethodName.Data(),"fCalculateQcumulants && !fPropagateErrorQC");
1838  }
1839  if(fCalculateQcumulants && fCalculateOnlySin)
1840  {
1841   Fatal(sMethodName.Data(),"fCalculateQcumulants && fCalculateOnlySin");
1842  }
1843  
1844  // e) Weights:
1845  for(Int_t rp=0;rp<2;rp++) // [RP,POI]
1846  {
1847   for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
1848   {
1849    if(fUseWeights[rp][ppe] && !fWeightsHist[rp][ppe])
1850    {
1851     Fatal(sMethodName.Data(),"fUseWeights[rp][ppe] && !fWeightsHist[rp][ppe], rp = %d, ppe = %d",rp,ppe);
1852    }
1853   }
1854  }
1855
1856  // f) Differential correlations:
1857  if(fCalculateDiffCorrelations && !fUseDefaultBinning && (fnDiffBins < 1 || !fRangesDiffBins))
1858  {
1859   Fatal(sMethodName.Data(),"fCalculateDiffCorrelations && !fUseDefaultBinning && (fnDiffBins < 1 || !fRangesDiffBins)"); 
1860  }
1861  if(fCalculateDiffCorrelations && !(fCalculateDiffCos || fCalculateDiffSin))
1862  {
1863   Fatal(sMethodName.Data(),"fCalculateDiffCorrelations && !(fCalculateDiffCos || fCalculateDiffSin)"); 
1864  }
1865
1866  // g) Nested loops:
1867  if(fCrossCheckDiffWithNestedLoops && (1 == fCrossCheckDiffCSCOBN[0] && !fCalculateDiffSin))
1868  {
1869   Fatal(sMethodName.Data(),"fCrossCheckDiffWithNestedLoops && (1 == fCrossCheckDiffCSCOBN[0] && !CalculateDiffSin)"); 
1870  }
1871  if(fCrossCheckDiffWithNestedLoops && (0 == fCrossCheckDiffCSCOBN[0] && !fCalculateDiffCos))
1872  {
1873   Fatal(sMethodName.Data(),"fCrossCheckDiffWithNestedLoops && (0 == fCrossCheckDiffCSCOBN[0] && !CalculateDiffCos)"); 
1874  }
1875
1876  // h) Dump the points:
1877  if(fDumpThePoints && !fFillMultDistributionsHist)
1878  {
1879   Fatal(sMethodName.Data(),"if(fDumpThePoints && !fFillMultDistributionsHist)"); 
1880  }
1881  if(fDumpThePoints && fMaxNoEventsPerFile <= 0)
1882  {
1883   Fatal(sMethodName.Data(),"if(fDumpThePoints && fMaxNoEventsPerFile <= 0)"); 
1884  }
1885
1886 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckSettings()
1887
1888 //=======================================================================================================================
1889
1890 void AliFlowAnalysisWithMultiparticleCorrelations::BookAndNestAllLists()
1891 {
1892  // Book and nest all lists nested in the base list fHistList.
1893
1894  // a) Book and nest lists for control histograms;
1895  // b) Book and nest lists for Q-vectors;
1896  // c) Book and nest lists for correlations;
1897  // d) Book and nest lists for e-b-e cumulants;
1898  // e) Book and nest lists for weights;
1899  // f) Book and nest lists for nested loops;
1900  // g) Book and nest lists for 'standard candles';
1901  // h) Book and nest lists for Q-cumulants;
1902  // i) Book and nest lists for differential correlations.
1903
1904  // a) Book and nest lists for control histograms:
1905  fControlHistogramsList = new TList();
1906  fControlHistogramsList->SetName("Control Histograms");
1907  fControlHistogramsList->SetOwner(kTRUE);
1908  fHistList->Add(fControlHistogramsList);
1909
1910  // b) Book and nest lists for Q-vectors:
1911  fQvectorList = new TList();
1912  fQvectorList->SetName("Q-vectors");
1913  fQvectorList->SetOwner(kTRUE);
1914  fHistList->Add(fQvectorList);
1915
1916  // c) Book and nest lists for correlations:
1917  fCorrelationsList = new TList();
1918  fCorrelationsList->SetName("Correlations");
1919  fCorrelationsList->SetOwner(kTRUE);
1920  fHistList->Add(fCorrelationsList);
1921
1922  // d) Book and nest lists for e-b-e cumulants:
1923  fEbECumulantsList = new TList();
1924  fEbECumulantsList->SetName("E-b-e Cumulants");
1925  fEbECumulantsList->SetOwner(kTRUE);
1926  fHistList->Add(fEbECumulantsList);
1927
1928  // e) Book and nest lists for weights:
1929  fWeightsList = new TList();
1930  fWeightsList->SetName("Weights");
1931  fWeightsList->SetOwner(kTRUE);
1932  fHistList->Add(fWeightsList);
1933
1934  // f) Book and nest lists for nested loops:
1935  fNestedLoopsList = new TList();
1936  fNestedLoopsList->SetName("Nested Loops");
1937  fNestedLoopsList->SetOwner(kTRUE);
1938  fHistList->Add(fNestedLoopsList);
1939
1940  // g) Book and nest lists for 'standard candles':
1941  fStandardCandlesList = new TList();
1942  fStandardCandlesList->SetName("Standard Candles");
1943  fStandardCandlesList->SetOwner(kTRUE);
1944  fHistList->Add(fStandardCandlesList);
1945
1946  // h) Book and nest lists for Q-cumulants:
1947  fQcumulantsList = new TList();
1948  fQcumulantsList->SetName("Q-cumulants");
1949  fQcumulantsList->SetOwner(kTRUE);
1950  fHistList->Add(fQcumulantsList);
1951
1952  // i) Book and nest lists for differential correlations:
1953  fDiffCorrelationsList = new TList();
1954  fDiffCorrelationsList->SetName("Differential Correlations");
1955  fDiffCorrelationsList->SetOwner(kTRUE);
1956  fHistList->Add(fDiffCorrelationsList);
1957
1958 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookAndNestAllLists()
1959
1960 //=======================================================================================================================
1961
1962 void AliFlowAnalysisWithMultiparticleCorrelations::WriteHistograms(TString outputFileName)
1963 {
1964  // Store the final results in output file <outputFileName>.root.
1965
1966  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
1967  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
1968
1969  delete output;
1970
1971 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::WriteHistograms(TString outputFileName)
1972
1973 //=======================================================================================================================
1974
1975 void AliFlowAnalysisWithMultiparticleCorrelations::WriteHistograms(TDirectoryFile *outputFileName)
1976 {
1977  // Store the final results in output file <outputFileName>.root.
1978
1979  outputFileName->Add(fHistList);
1980  outputFileName->Write(outputFileName->GetName(),TObject::kSingleKey);
1981
1982 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::WriteHistograms(TDirectoryFile *outputFileName)
1983
1984 //=======================================================================================================================
1985
1986 void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForControlHistograms()
1987 {
1988  // Book all the stuff for control histograms.
1989
1990  // a) Book the profile holding all the flags for control histograms;
1991  // b) Book all control histograms;
1992  //  b0) Book TH1D *fKinematicsHist[2][3]; 
1993  //  b1) Book TH1D *fMultDistributionsHist[3];
1994  //  b2) Book TH2D *fMultCorrelationsHist[3].
1995
1996  // a) Book the profile holding all the flags for control histograms: TBI stil incomplete 
1997  fControlHistogramsFlagsPro = new TProfile("fControlHistogramsFlagsPro","Flags and settings for control histograms",4,0,4);
1998  fControlHistogramsFlagsPro->SetTickLength(-0.01,"Y");
1999  fControlHistogramsFlagsPro->SetMarkerStyle(25);
2000  fControlHistogramsFlagsPro->SetLabelSize(0.04);
2001  fControlHistogramsFlagsPro->SetLabelOffset(0.02,"Y");
2002  fControlHistogramsFlagsPro->SetStats(kFALSE);
2003  fControlHistogramsFlagsPro->SetFillColor(kGray);
2004  fControlHistogramsFlagsPro->SetLineColor(kBlack);
2005  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistograms"); fControlHistogramsFlagsPro->Fill(0.5,fFillControlHistograms);
2006  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(2,"fFillKinematicsHist"); fControlHistogramsFlagsPro->Fill(1.5,fFillKinematicsHist);
2007  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(3,"fFillMultDistributionsHist"); fControlHistogramsFlagsPro->Fill(2.5,fFillMultDistributionsHist);
2008  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(4,"fFillMultCorrelationsHist"); fControlHistogramsFlagsPro->Fill(3.5,fFillMultCorrelationsHist);
2009  fControlHistogramsList->Add(fControlHistogramsFlagsPro);
2010
2011  if(!fFillControlHistograms){return;} // TBI is this safe? Well, perhaps it is if I can't implement it better...
2012
2013  // b) Book all control histograms: // TBI add setters for all these values
2014  //  b0) Book TH1D *fKinematicsHist[2][3]:
2015  TString name[2][3] = {{"RP,phi","RP,pt","RP,eta"},{"POI,phi","POI,pt","POI,eta"}}; // [RP,POI][phi,pt,eta]
2016  TString title[2] = {"Reference particles (RP)","Particles of interest (POI)"}; // [RP,POI]
2017  Int_t lineColor[2] = {kBlue,kRed}; // [RP,POI]
2018  Int_t fillColor[2] = {kBlue-10,kRed-10}; // [RP,POI]
2019  TString xAxisTitle[3] = {"#phi","p_{T}","#eta"}; // [phi,pt,eta]
2020  if(fFillKinematicsHist)
2021  {
2022   for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2023   {
2024    for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2025    {
2026     fKinematicsHist[rp][ppe] = new TH1D(name[rp][ppe].Data(),title[rp].Data(),fnBins[rp][ppe],fMin[rp][ppe],fMax[rp][ppe]);
2027     fKinematicsHist[rp][ppe]->GetXaxis()->SetTitle(xAxisTitle[ppe].Data());
2028     fKinematicsHist[rp][ppe]->SetLineColor(lineColor[rp]);
2029     fKinematicsHist[rp][ppe]->SetFillColor(fillColor[rp]);
2030     fKinematicsHist[rp][ppe]->SetMinimum(0.); 
2031     fControlHistogramsList->Add(fKinematicsHist[rp][ppe]);
2032    }
2033   }
2034  } // if(fFillKinematicsHist)
2035
2036  //  b1) Book TH1D *fMultDistributionsHist[3]: // TBI add setters for all these values
2037  TString nameMult[3] = {"Multiplicity (RP)","Multiplicity (POI)","Multiplicity (REF)"}; // [RP,POI,reference multiplicity]
2038  TString titleMult[3] = {"Reference particles (RP)","Particles of interest (POI)",""}; // [RP,POI,reference multiplicity]
2039  Int_t lineColorMult[3] = {kBlue,kRed,kGreen+2}; // [RP,POI,reference multiplicity]
2040  Int_t fillColorMult[3] = {kBlue-10,kRed-10,kGreen-10}; // [RP,POI,reference multiplicity]
2041  TString xAxisTitleMult[3] = {"Multiplicity (RP)","Multiplicity (POI)","Multiplicity (REF)"}; // [phi,pt,eta]
2042  if(fFillMultDistributionsHist)
2043  {
2044   for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2045   {
2046    fMultDistributionsHist[rprm] = new TH1D(nameMult[rprm].Data(),titleMult[rprm].Data(),fnBinsMult[rprm],fMinMult[rprm],fMaxMult[rprm]);
2047    fMultDistributionsHist[rprm]->GetXaxis()->SetTitle(xAxisTitleMult[rprm].Data());
2048    fMultDistributionsHist[rprm]->SetLineColor(lineColorMult[rprm]);
2049    fMultDistributionsHist[rprm]->SetFillColor(fillColorMult[rprm]);
2050    fControlHistogramsList->Add(fMultDistributionsHist[rprm]);
2051   } // for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2052  } // if(fFillMultDistributionsHist)
2053
2054  //  b2) Book TH2I *fMultCorrelationsHist[3]: 
2055  if(fFillMultCorrelationsHist)
2056  {
2057   // ...
2058   fMultCorrelationsHist[0] = new TH2I("Multiplicity (RP vs. POI)","Multiplicity (RP vs. POI)",fnBinsMult[0],fMinMult[0],fMaxMult[0],fnBinsMult[1],fMinMult[1],fMaxMult[1]);
2059   fMultCorrelationsHist[0]->GetXaxis()->SetTitle(xAxisTitleMult[0].Data());
2060   fMultCorrelationsHist[0]->GetYaxis()->SetTitle(xAxisTitleMult[1].Data());
2061   fControlHistogramsList->Add(fMultCorrelationsHist[0]);
2062   // ...
2063   fMultCorrelationsHist[1] = new TH2I("Multiplicity (RP vs. REF)","Multiplicity (RP vs. REF)",fnBinsMult[0],fMinMult[0],fMaxMult[0],fnBinsMult[2],fMinMult[2],fMaxMult[2]);
2064   fMultCorrelationsHist[1]->GetXaxis()->SetTitle(xAxisTitleMult[0].Data());
2065   fMultCorrelationsHist[1]->GetYaxis()->SetTitle(xAxisTitleMult[2].Data());
2066   fControlHistogramsList->Add(fMultCorrelationsHist[1]);
2067   // ...
2068   fMultCorrelationsHist[2] = new TH2I("Multiplicity (POI vs. REF)","Multiplicity (POI vs. REF)",fnBinsMult[1],fMinMult[1],fMaxMult[1],fnBinsMult[2],fMinMult[2],fMaxMult[2]);
2069   fMultCorrelationsHist[2]->GetXaxis()->SetTitle(xAxisTitleMult[1].Data());
2070   fMultCorrelationsHist[2]->GetYaxis()->SetTitle(xAxisTitleMult[2].Data());
2071   fControlHistogramsList->Add(fMultCorrelationsHist[2]);
2072  } // if(fFillMultCorrelationsHist){
2073
2074 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForControlHistograms()
2075
2076 //=======================================================================================================================
2077
2078 void AliFlowAnalysisWithMultiparticleCorrelations::FillControlHistograms(AliFlowEventSimple *anEvent)
2079 {
2080  // Fill control histograms. 
2081  // a) Fill TH1D *fKinematicsHist[2][3];
2082  // b) Fill TH1D *fMultDistributionsHist[3]; 
2083  // c) Fill TH2D *fMultCorrelationsHist[3].  
2084
2085  // a) Fill TH1D *fKinematicsHist[2][3]:
2086  if(fFillKinematicsHist)
2087  {
2088   Int_t nTracks = anEvent->NumberOfTracks(); // TBI shall I promote this to data member?
2089   for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2090   {
2091    AliFlowTrackSimple *pTrack = anEvent->GetTrack(t);
2092    if(!pTrack){printf("\n AAAARGH: pTrack is NULL in MPC::FCH() !!!!");continue;}
2093    if(pTrack)
2094    {
2095     Double_t dPhi = pTrack->Phi(); 
2096     //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
2097     //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
2098     Double_t dPt = pTrack->Pt();
2099     Double_t dEta = pTrack->Eta();
2100     Double_t dPhiPtEta[3] = {dPhi,dPt,dEta};
2101     for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2102     {
2103      for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2104      {
2105       if((0==rp && pTrack->InRPSelection()) || (1==rp && pTrack->InPOISelection())) // TBI 
2106       { 
2107        fKinematicsHist[rp][ppe]->Fill(dPhiPtEta[ppe]);
2108       }
2109      } // for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2110     } // for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2111    } // if(pTrack)  
2112   } // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2113  } // if(fFillKinematicsHist)
2114
2115  // b) Fill TH1D *fMultDistributionsHist[3]: 
2116  Double_t dMultRP = anEvent->GetNumberOfRPs(); // TBI shall I promote these 3 variables into data members? 
2117  Double_t dMultPOI = anEvent->GetNumberOfPOIs();
2118  Double_t dMultREF = anEvent->GetReferenceMultiplicity();
2119  Double_t dMult[3] = {dMultRP,dMultPOI,dMultREF};
2120  for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2121  {
2122   if(fFillMultDistributionsHist){fMultDistributionsHist[rprm]->Fill(dMult[rprm]);}      
2123  } 
2124
2125  // c) Fill TH2I *fMultCorrelationsHist[3]:  
2126  if(fFillMultCorrelationsHist)
2127  {
2128   fMultCorrelationsHist[0]->Fill((Int_t)dMultRP,(Int_t)dMultPOI); // RP vs. POI
2129   fMultCorrelationsHist[1]->Fill((Int_t)dMultRP,(Int_t)dMultREF); // RP vs. refMult
2130   fMultCorrelationsHist[2]->Fill((Int_t)dMultPOI,(Int_t)dMultREF); // POI vs. refMult
2131  } // if(fFillMultCorrelationsHist)
2132
2133 } // void AliFlowAnalysisWithMultiparticleCorrelations::FillControlHistograms(AliFlowEventSimple *anEvent)
2134
2135 //=======================================================================================================================
2136
2137 void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForControlHistograms()
2138 {
2139  // Initialize all arrays for control histograms.
2140
2141  // a) Initialize TH1D *fKinematicsHist[2][3];
2142  // b) Initialize TH1D *fMultDistributionsHist[3]; 
2143  // c) Initialize TH2D *fMultCorrelationsHist[3];  
2144  // d) Initialize default binning values for fKinematicsHist[2][3];
2145  // e) Initialize default binning values for fMultCorrelationsHist[3].
2146  
2147  // a) Initialize TH1D *fKinematicsHist[2][3]:
2148  for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2149  {
2150   for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2151   {
2152    fKinematicsHist[rp][ppe] = NULL;
2153   } 
2154  } 
2155
2156  // b) Initialize TH1D *fMultDistributionsHist[3]:
2157  for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2158  {
2159   fMultDistributionsHist[rprm] = NULL;      
2160  } 
2161
2162  // c) Initialize TH2I *fMultCorrelationsHist[3]: 
2163  for(Int_t r=0;r<3;r++) // [RP vs. POI, RP vs. refMult, POI vs. refMult]  
2164  {
2165   fMultCorrelationsHist[r] = NULL; 
2166  }
2167
2168  // d) Initialize default binning values for fKinematicsHist[2][3]:
2169  // nBins:
2170  fnBins[0][0] = 360;  // [RP][phi]
2171  fnBins[0][1] = 1000; // [RP][pt]
2172  fnBins[0][2] = 1000; // [RP][eta]
2173  fnBins[1][0] = 360;  // [POI][phi]
2174  fnBins[1][1] = 1000; // [POI][pt]
2175  fnBins[1][2] = 1000; // [POI][eta]
2176  // Min:
2177  fMin[0][0] = 0.;  // [RP][phi]
2178  fMin[0][1] = 0.;  // [RP][pt]
2179  fMin[0][2] = -1.; // [RP][eta]
2180  fMin[1][0] = 0.;  // [POI][phi]
2181  fMin[1][1] = 0.;  // [POI][pt]
2182  fMin[1][2] = -1.; // [POI][eta]
2183  // Max:
2184  fMax[0][0] = TMath::TwoPi(); // [RP][phi]
2185  fMax[0][1] = 10.;            // [RP][pt]
2186  fMax[0][2] = 1.;             // [RP][eta]
2187  fMax[1][0] = TMath::TwoPi(); // [POI][phi]
2188  fMax[1][1] = 10.;            // [POI][pt]
2189  fMax[1][2] = 1.;             // [POI][eta]
2190
2191  // e) Initialize default binning values for fMultCorrelationsHist[3]:
2192  // nBins:
2193  fnBinsMult[0] = 3000; // [RP]
2194  fnBinsMult[1] = 3000; // [POI]
2195  fnBinsMult[2] = 3000; // [REF]
2196  // Min:
2197  fMinMult[0] = 0.; // [RP]
2198  fMinMult[1] = 0.; // [POI]
2199  fMinMult[2] = 0.; // [REF]
2200  // Max:
2201  fMaxMult[0] = 3000.; // [RP]
2202  fMaxMult[1] = 3000.; // [POI]
2203  fMaxMult[2] = 3000.; // [REF]
2204
2205 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForControlHistograms()
2206
2207 //=======================================================================================================================
2208
2209 void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForQvector()
2210 {
2211  // Book all the stuff for Q-vector.
2212
2213  // a) Book the profile holding all the flags for Q-vector;
2214  // ...
2215
2216  // a) Book the profile holding all the flags for Q-vector:
2217  fQvectorFlagsPro = new TProfile("fQvectorFlagsPro","Flags for Q-vectors",2,0,2);
2218  fQvectorFlagsPro->SetTickLength(-0.01,"Y");
2219  fQvectorFlagsPro->SetMarkerStyle(25);
2220  fQvectorFlagsPro->SetLabelSize(0.03);
2221  fQvectorFlagsPro->SetLabelOffset(0.02,"Y");
2222  fQvectorFlagsPro->SetStats(kFALSE);
2223  fQvectorFlagsPro->SetFillColor(kGray);
2224  fQvectorFlagsPro->SetLineColor(kBlack);
2225  fQvectorFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateQvector"); fQvectorFlagsPro->Fill(0.5,fCalculateQvector); 
2226  fQvectorFlagsPro->GetXaxis()->SetBinLabel(2,"fCalculateDiffQvectors"); fQvectorFlagsPro->Fill(1.5,fCalculateDiffQvectors); 
2227  fQvectorList->Add(fQvectorFlagsPro);
2228
2229  // ...
2230
2231 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForQvector()
2232
2233 //=======================================================================================================================
2234
2235 void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForCorrelations()
2236 {
2237  // Book all the stuff for correlations.
2238
2239  // TBI this method can be implemented in a much more civilised way. 
2240
2241  // a) Book the profile holding all the flags for correlations;
2242  // b) Book TProfile *fCorrelationsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]).
2243
2244  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForCorrelations()";
2245
2246  // a) Book the profile holding all the flags for correlations:
2247  fCorrelationsFlagsPro = new TProfile("fCorrelationsFlagsPro","Flags for correlations",13,0,13);
2248  fCorrelationsFlagsPro->SetTickLength(-0.01,"Y");
2249  fCorrelationsFlagsPro->SetMarkerStyle(25);
2250  fCorrelationsFlagsPro->SetLabelSize(0.03);
2251  fCorrelationsFlagsPro->SetLabelOffset(0.02,"Y");
2252  fCorrelationsFlagsPro->SetStats(kFALSE);
2253  fCorrelationsFlagsPro->SetFillColor(kGray);
2254  fCorrelationsFlagsPro->SetLineColor(kBlack);
2255  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateCorrelations"); fCorrelationsFlagsPro->Fill(0.5,fCalculateCorrelations); 
2256  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(2,"fMaxHarmonic"); fCorrelationsFlagsPro->Fill(1.5,fMaxHarmonic); 
2257  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(3,"fMaxCorrelator"); fCorrelationsFlagsPro->Fill(2.5,fMaxCorrelator); 
2258  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(4,"fCalculateIsotropic"); fCorrelationsFlagsPro->Fill(3.5,fCalculateIsotropic); 
2259  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(5,"fCalculateSame"); fCorrelationsFlagsPro->Fill(4.5,fCalculateSame); 
2260  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(6,"fSkipZeroHarmonics"); fCorrelationsFlagsPro->Fill(5.5,fSkipZeroHarmonics); 
2261  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(7,"fCalculateSameIsotropic"); fCorrelationsFlagsPro->Fill(6.5,fCalculateSameIsotropic); 
2262  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(8,"fCalculateAll"); fCorrelationsFlagsPro->Fill(7.5,fCalculateAll); 
2263  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(9,"fDontGoBeyond"); fCorrelationsFlagsPro->Fill(8.5,fDontGoBeyond); 
2264  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(10,"fCalculateOnlyForHarmonicQC"); fCorrelationsFlagsPro->Fill(9.5,fCalculateOnlyForHarmonicQC); 
2265  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(11,"fCalculateOnlyForSC"); fCorrelationsFlagsPro->Fill(10.5,fCalculateOnlyForSC); 
2266  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(12,"fCalculateOnlyCos"); fCorrelationsFlagsPro->Fill(11.5,fCalculateOnlyCos); 
2267  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(13,"fCalculateOnlySin"); fCorrelationsFlagsPro->Fill(12.5,fCalculateOnlySin);
2268  fCorrelationsList->Add(fCorrelationsFlagsPro);
2269
2270  if(!fCalculateCorrelations){return;} // TBI is this safe enough? 
2271
2272  // b) Book TProfile *fCorrelationsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]): // TBI hardwired 8, shall be fMaxCorrelator
2273  cout<<" => Booking TProfile *fCorrelationsPro[2][8]..."<<endl;
2274  TString sCosSin[2] = {"Cos","Sin"};
2275  Int_t markerColor[2] = {kBlue,kRed};
2276  Int_t markerStyle[2] = {kFullSquare,kFullSquare};
2277  Int_t nBins[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
2278  Int_t nBinsTitle[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
2279  Int_t nToBeFilled[8] = {0,0,0,0,0,0,0,0}; // TBI hardwired 8, shall be fMaxCorrelator
2280  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2281  {
2282   // Implementing \binom{n+k-1}{k}, which is the resulting number of sets obtained
2283   // after sampling n starting elements into k subsets, repetitions allowed.
2284   // In my case, n=2*fMaxHarmonic+1, k=c+1, hence:
2285   nBins[c] = (Int_t)(TMath::Factorial(2*fMaxHarmonic+1+c+1-1)
2286            / (TMath::Factorial(2*fMaxHarmonic+1-1)*TMath::Factorial(c+1)));
2287   nBinsTitle[c] = nBins[c];
2288   if(c>=fDontGoBeyond){nBins[c]=1;} // TBI is this really safe? 
2289  } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
2290  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2291  {
2292   if(fCalculateOnlyCos && 1==cs){continue;}
2293   else if(fCalculateOnlySin && 0==cs){continue;}
2294   for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2295   {
2296    fCorrelationsPro[cs][c] = new TProfile(Form("%dpCorrelations%s",c+1,sCosSin[cs].Data()),"",nBins[c],0.,1.*nBins[c]);
2297    fCorrelationsPro[cs][c]->Sumw2();
2298    fCorrelationsPro[cs][c]->SetStats(kFALSE);
2299    fCorrelationsPro[cs][c]->SetMarkerColor(markerColor[cs]);
2300    fCorrelationsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
2301    fCorrelationsList->Add(fCorrelationsPro[cs][c]);
2302   } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
2303  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2304  // Set all bin labels: TBI this can be implemented better, most likely...
2305  Int_t binNo[2][8]; 
2306  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2307  {
2308   if(fCalculateOnlyCos && 1==cs){continue;}
2309   else if(fCalculateOnlySin && 0==cs){continue;}
2310   for(Int_t c=0;c<fMaxCorrelator;c++)
2311   {
2312    binNo[cs][c] = 1;
2313   } 
2314  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2315  
2316  for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++) 
2317  {
2318   cout<< Form("    Patience, this takes some time... n1 = %d/%d\r",n1+fMaxHarmonic,2*fMaxHarmonic)<<flush; // TBI
2319   if(fSkipZeroHarmonics && 0==n1){continue;}
2320   if(fCalculateOnlyForHarmonicQC && TMath::Abs(n1) != fHarmonicQC){continue;}
2321   if(fCalculateAll)
2322   {
2323    for(Int_t cs=0;cs<2;cs++) 
2324    {
2325     if(fCalculateOnlyCos && 1==cs){continue;}
2326     else if(fCalculateOnlySin && 0==cs){continue;}
2327     fCorrelationsPro[cs][0]->GetXaxis()->SetBinLabel(binNo[cs][0]++,Form("%s(%d)",sCosSin[cs].Data(),n1));
2328    } // for(Int_t cs=0;cs<2;cs++) 
2329    nToBeFilled[0]++; 
2330   }
2331   if(1==fDontGoBeyond){continue;}
2332   for(Int_t n2=n1;n2<=fMaxHarmonic;n2++) 
2333   {
2334    if(fSkipZeroHarmonics && 0==n2){continue;}
2335    if(fCalculateOnlyForHarmonicQC && TMath::Abs(n2) != fHarmonicQC){continue;}
2336    if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2)) 
2337       || (fCalculateSameIsotropic && 0==n1+n2 && TMath::Abs(n1)==TMath::Abs(n2)) 
2338       || (fCalculateOnlyForHarmonicQC && 0==n1+n2)
2339       || (fCalculateOnlyForSC && 0==n1+n2))
2340    {  
2341     for(Int_t cs=0;cs<2;cs++) 
2342     {
2343      if(fCalculateOnlyCos && 1==cs){continue;}
2344      else if(fCalculateOnlySin && 0==cs){continue;}
2345      fCorrelationsPro[cs][1]->GetXaxis()->SetBinLabel(binNo[cs][1]++,Form("%s(%d,%d)",sCosSin[cs].Data(),n1,n2));
2346     } // for(Int_t cs=0;cs<2;cs++) 
2347     nToBeFilled[1]++; 
2348    }
2349    if(2==fDontGoBeyond){continue;}
2350    for(Int_t n3=n2;n3<=fMaxHarmonic;n3++) 
2351    {
2352     if(fSkipZeroHarmonics && 0==n3){continue;}
2353     if(fCalculateOnlyForHarmonicQC && TMath::Abs(n3) != fHarmonicQC){continue;}
2354     if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3))
2355        || (fCalculateSameIsotropic && 0==n1+n2+n3 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)) 
2356        || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3))
2357     {  
2358      for(Int_t cs=0;cs<2;cs++) 
2359      {
2360       if(fCalculateOnlyCos && 1==cs){continue;}
2361       else if(fCalculateOnlySin && 0==cs){continue;}
2362       fCorrelationsPro[cs][2]->GetXaxis()->SetBinLabel(binNo[cs][2]++,Form("%s(%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3));
2363      } // for(Int_t cs=0;cs<2;cs++) 
2364      nToBeFilled[2]++; 
2365     }
2366     if(3==fDontGoBeyond){continue;}
2367     for(Int_t n4=n3;n4<=fMaxHarmonic;n4++) 
2368     {
2369      if(fSkipZeroHarmonics && 0==n4){continue;}
2370      if(fCalculateOnlyForHarmonicQC && TMath::Abs(n4) != fHarmonicQC){continue;}
2371      if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4))
2372        || (fCalculateSameIsotropic && 0==n1+n2+n3+n4 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)) 
2373        || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4)
2374        || (fCalculateOnlyForSC && (0==n1+n4 && 0==n2+n3 && n1 != n2 && n3 != n4)))
2375      {   
2376       for(Int_t cs=0;cs<2;cs++) 
2377       {
2378        if(fCalculateOnlyCos && 1==cs){continue;}
2379        else if(fCalculateOnlySin && 0==cs){continue;}
2380        fCorrelationsPro[cs][3]->GetXaxis()->SetBinLabel(binNo[cs][3]++,Form("%s(%d,%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3,n4));
2381       } // for(Int_t cs=0;cs<2;cs++) 
2382       nToBeFilled[3]++; 
2383      } 
2384      if(4==fDontGoBeyond){continue;}
2385      for(Int_t n5=n4;n5<=fMaxHarmonic;n5++) 
2386      {
2387       if(fSkipZeroHarmonics && 0==n5){continue;}
2388       if(fCalculateOnlyForHarmonicQC && TMath::Abs(n5) != fHarmonicQC){continue;}
2389       if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5) 
2390          || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5))
2391          || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5)) 
2392          || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5))
2393       {   
2394        for(Int_t cs=0;cs<2;cs++) 
2395        {
2396         if(fCalculateOnlyCos && 1==cs){continue;}
2397         else if(fCalculateOnlySin && 0==cs){continue;}
2398         fCorrelationsPro[cs][4]->GetXaxis()->SetBinLabel(binNo[cs][4]++,Form("%s(%d,%d,%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3,n4,n5));
2399        } // for(Int_t cs=0;cs<2;cs++) 
2400        nToBeFilled[4]++; 
2401       }
2402       if(5==fDontGoBeyond){continue;}
2403       for(Int_t n6=n5;n6<=fMaxHarmonic;n6++) 
2404       {
2405        if(fSkipZeroHarmonics && 0==n6){continue;}
2406        if(fCalculateOnlyForHarmonicQC && TMath::Abs(n6) != fHarmonicQC){continue;}
2407        if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6)  
2408           || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2409               && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6))
2410           || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
2411               && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6)) 
2412           || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5+n6))
2413        {   
2414         for(Int_t cs=0;cs<2;cs++) 
2415         {
2416          if(fCalculateOnlyCos && 1==cs){continue;}
2417          else if(fCalculateOnlySin && 0==cs){continue;}
2418          fCorrelationsPro[cs][5]->GetXaxis()->SetBinLabel(binNo[cs][5]++,Form("%s(%d,%d,%d,%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3,n4,n5,n6));         
2419         } // for(Int_t cs=0;cs<2;cs++) 
2420         nToBeFilled[5]++; 
2421        }
2422        if(6==fDontGoBeyond){continue;}
2423        for(Int_t n7=n6;n7<=fMaxHarmonic;n7++) 
2424        {
2425         if(fSkipZeroHarmonics && 0==n7){continue;}
2426         if(fCalculateOnlyForHarmonicQC && TMath::Abs(n7) != fHarmonicQC){continue;}
2427         if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7) 
2428            || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2429                && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7))
2430            || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2431                && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7)) 
2432            || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5+n6+n7))
2433         {   
2434          for(Int_t cs=0;cs<2;cs++) 
2435          {
2436           if(fCalculateOnlyCos && 1==cs){continue;}
2437           else if(fCalculateOnlySin && 0==cs){continue;}
2438           fCorrelationsPro[cs][6]->GetXaxis()->SetBinLabel(binNo[cs][6]++,Form("%s(%d,%d,%d,%d,%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3,n4,n5,n6,n7));
2439          } // for(Int_t cs=0;cs<2;cs++) 
2440          nToBeFilled[6]++; 
2441         }
2442         if(7==fDontGoBeyond){continue;}
2443         for(Int_t n8=n7;n8<=fMaxHarmonic;n8++) 
2444         {
2445          if(fSkipZeroHarmonics && 0==n8){continue;}
2446          if(fCalculateOnlyForHarmonicQC && TMath::Abs(n8) != fHarmonicQC){continue;}
2447          if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8) 
2448             || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2449                 && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7) && TMath::Abs(n1)==TMath::Abs(n8))
2450             || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
2451                 && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7) 
2452                 && TMath::Abs(n1)==TMath::Abs(n8))
2453             || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5+n6+n7+n8))
2454          {    
2455           for(Int_t cs=0;cs<2;cs++) 
2456           {
2457            if(fCalculateOnlyCos && 1==cs){continue;}
2458            else if(fCalculateOnlySin && 0==cs){continue;}
2459            fCorrelationsPro[cs][7]->GetXaxis()->SetBinLabel(binNo[cs][7]++,Form("%s(%d,%d,%d,%d,%d,%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3,n4,n5,n6,n7,n8));
2460           } // for(Int_t cs=0;cs<2;cs++) 
2461           nToBeFilled[7]++; 
2462          }
2463         } // for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
2464        } // for(Int_t n7=n6;n7<=fMaxHarmonic;n7++) 
2465       } // for(Int_t n6=n5;n6<=fMaxHarmonic;n6++) 
2466      } // for(Int_t n5=n4;n5<=fMaxHarmonic;n5++) 
2467     } // for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)   
2468    } // for(Int_t n3=n2;n3<=fMaxHarmonic;n3++) 
2469   } // for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
2470  } // for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++) 
2471
2472  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2473  {
2474   if(fCalculateOnlyCos && 1==cs){continue;}
2475   else if(fCalculateOnlySin && 0==cs){continue;}
2476   for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2477   {
2478    fCorrelationsPro[cs][c]->SetTitle(Form("%d-p correlations, %s terms, %d/%d in total",c+1,sCosSin[cs].Data(),nToBeFilled[c],nBinsTitle[c]));
2479    fCorrelationsPro[cs][c]->GetXaxis()->SetRangeUser(0.,fCorrelationsPro[cs][c]->GetBinLowEdge(nToBeFilled[c]+1));
2480   }
2481  } 
2482  cout<<"    Booked.                                           "<<endl; // TBI 
2483
2484 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForCorrelations()
2485
2486 //=======================================================================================================================
2487
2488 void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForDiffCorrelations()
2489 {
2490  // Book all the stuff for differential correlations.
2491
2492  // a) Book the profile holding all the flags for differential correlations;
2493  // b) Book TProfile *fDiffCorrelationsPro[2][4] ([0=cos,1=sin][1p,2p,3p,4p]).
2494
2495  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForDiffCorrelations()";
2496
2497  // a) Book the profile holding all the flags for differential correlations:
2498  fDiffCorrelationsFlagsPro = new TProfile("fDiffCorrelationsFlagsPro","Flags for differential correlations",5,0,5);
2499  fDiffCorrelationsFlagsPro->SetTickLength(-0.01,"Y");
2500  fDiffCorrelationsFlagsPro->SetMarkerStyle(25);
2501  fDiffCorrelationsFlagsPro->SetLabelSize(0.03);
2502  fDiffCorrelationsFlagsPro->SetLabelOffset(0.02,"Y");
2503  fDiffCorrelationsFlagsPro->SetStats(kFALSE);
2504  fDiffCorrelationsFlagsPro->SetFillColor(kGray);
2505  fDiffCorrelationsFlagsPro->SetLineColor(kBlack);
2506  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateDiffCorrelations"); fDiffCorrelationsFlagsPro->Fill(0.5,fCalculateDiffCorrelations); 
2507  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(2,"fCalculateDiffCos"); fDiffCorrelationsFlagsPro->Fill(1.5,fCalculateDiffCos); 
2508  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(3,"fCalculateDiffSin"); fDiffCorrelationsFlagsPro->Fill(2.5,fCalculateDiffSin); 
2509  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(4,"fCalculateDiffCorrelationsVsPt"); fDiffCorrelationsFlagsPro->Fill(3.5,fCalculateDiffCorrelationsVsPt); 
2510  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(5,"fUseDefaultBinning"); fDiffCorrelationsFlagsPro->Fill(4.5,fUseDefaultBinning); 
2511  fDiffCorrelationsList->Add(fDiffCorrelationsFlagsPro);
2512
2513  // b) Book TProfile *fDiffCorrelationsPro[2][4] ([0=cos,1=sin][1p,2p,3p,4p]):
2514  Bool_t fDiffStore[2][4] = {{0,1,1,1},{0,0,0,0}}; // store or not TBI promote to data member, and implement setter perhaps  
2515  Int_t markerColor[2] = {kRed,kGreen};
2516  Int_t markerStyle[2] = {kFullSquare,kOpenSquare};
2517  TString sCosSin[2] = {"Cos","Sin"};
2518  TString sLabel[4] = {Form("%d",fDiffHarmonics[0][0]),
2519                       Form("%d,%d",fDiffHarmonics[1][0],fDiffHarmonics[1][1]),
2520                       Form("%d,%d,%d",fDiffHarmonics[2][0],fDiffHarmonics[2][1],fDiffHarmonics[2][2]),
2521                       Form("%d,%d,%d,%d",fDiffHarmonics[3][0],fDiffHarmonics[3][1],fDiffHarmonics[3][2],fDiffHarmonics[3][3])};
2522
2523  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2524  {
2525   if(!fCalculateDiffCos && 0==cs){continue;}
2526   if(!fCalculateDiffSin && 1==cs){continue;}
2527
2528   for(Int_t c=0;c<4;c++) // [1p,2p,3p,4p]
2529   {
2530    if(fCalculateDiffCorrelationsVsPt)
2531    {
2532     if(fUseDefaultBinning)
2533     {
2534      // vs pt, default binning:  
2535      fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"pt"),
2536                                                 Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
2537                                                 100,0.,10.);
2538     } else // if(fUseDefaultBinning)
2539       {
2540        // vs pt, non-default binning:
2541        fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"pt"),
2542                                                   Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
2543                                                   fnDiffBins,fRangesDiffBins);
2544       }// else // if(fUseDefaultBinning) 
2545       fDiffCorrelationsPro[cs][c]->Sumw2();
2546       fDiffCorrelationsPro[cs][c]->SetStats(kFALSE);
2547       fDiffCorrelationsPro[cs][c]->SetMarkerColor(markerColor[cs]);
2548       fDiffCorrelationsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
2549       fDiffCorrelationsPro[cs][c]->GetXaxis()->SetTitle("p_{T}");
2550       if(fDiffStore[cs][c]){fDiffCorrelationsList->Add(fDiffCorrelationsPro[cs][c]);}
2551    } else // if(fCalculateDiffCorrelationsVsPt)
2552      {
2553       if(fUseDefaultBinning)
2554       {
2555        // vs eta, default binning:
2556        fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"eta"),
2557                                                   Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
2558                                                   100,-1.,1.);
2559       } else // if(fUseDefaultBinning)
2560         {
2561          // vs eta, non-default binning:
2562          fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"eta"),
2563                                                     Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
2564                                                     fnDiffBins,fRangesDiffBins);
2565         } // else // if(fUseDefaultBinning)
2566         fDiffCorrelationsPro[cs][c]->Sumw2();
2567         fDiffCorrelationsPro[cs][c]->SetStats(kFALSE);
2568         fDiffCorrelationsPro[cs][c]->SetMarkerColor(markerColor[cs]);
2569         fDiffCorrelationsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
2570         fDiffCorrelationsPro[cs][c]->GetXaxis()->SetTitle("#eta");
2571         if(fDiffStore[cs][c]){fDiffCorrelationsList->Add(fDiffCorrelationsPro[cs][c]);}
2572      } // else // if(fCalculateDiffCorrelationsVsPt)
2573   } // for(Int_t c=0;c<4;c++) // [1p,2p,3p,4p]
2574  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2575
2576 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForDiffCorrelations()
2577
2578 //=======================================================================================================================
2579
2580 void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForEbECumulants()
2581 {
2582  // Book all the stuff for event-by-event cumulants.
2583
2584  // a) Book the profile holding all the flags for e-b-e cumulants;
2585  // b) Book TProfile *fEbECumulantsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]).
2586
2587  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForEbECumulants()";
2588
2589  // a) Book the profile holding all the flags for e-b-e cumulants:
2590  fEbECumulantsFlagsPro = new TProfile("fEbECumulantsFlagsPro","Flags for e-b-e cumulants",1,0,1);
2591  fEbECumulantsFlagsPro->SetTickLength(-0.01,"Y");
2592  fEbECumulantsFlagsPro->SetMarkerStyle(25);
2593  fEbECumulantsFlagsPro->SetLabelSize(0.03);
2594  fEbECumulantsFlagsPro->SetLabelOffset(0.02,"Y");
2595  fEbECumulantsFlagsPro->SetStats(kFALSE);
2596  fEbECumulantsFlagsPro->SetFillColor(kGray);
2597  fEbECumulantsFlagsPro->SetLineColor(kBlack);
2598  fEbECumulantsFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateEbECumulants"); fEbECumulantsFlagsPro->Fill(0.5,fCalculateEbECumulants); 
2599  fEbECumulantsList->Add(fEbECumulantsFlagsPro);
2600
2601  if(!fCalculateEbECumulants){return;} // TBI is this safe enough? 
2602
2603  // b) Book TProfile *fEbECumulantsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]):
2604  TString sCosSin[2] = {"Cos","Sin"};
2605  Int_t markerColor[2] = {kBlue,kRed};
2606  Int_t markerStyle[2] = {kFullSquare,kFullSquare};
2607  Int_t nBins[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
2608  Int_t nBinsTitle[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
2609  Int_t nToBeFilled[8] = {0,0,0,0,0,0,0,0}; // TBI hardwired 8, shall be fMaxCorrelator
2610  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2611  {
2612   // Implementing \binom{n+k-1}{k}, which is the resulting number of sets obtained
2613   // after sampling n starting elements into k subsets, repetitions allowed.
2614   // In my case, n=2*fMaxHarmonic+1, k=c+1, hence:
2615   nBins[c] = (Int_t)(TMath::Factorial(2*fMaxHarmonic+1+c+1-1)
2616            / (TMath::Factorial(2*fMaxHarmonic+1-1)*TMath::Factorial(c+1)));
2617   nBinsTitle[c] = nBins[c];
2618   if(c>=fDontGoBeyond){nBins[c]=1;} // TBI a bit of spaghetti here...
2619  } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
2620  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2621  {
2622   for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2623   {
2624    fEbECumulantsPro[cs][c] = new TProfile(Form("%dpEbECumulants%s",c+1,sCosSin[cs].Data()),"",nBins[c],0.,1.*nBins[c]);
2625    fEbECumulantsPro[cs][c]->Sumw2();
2626    fEbECumulantsPro[cs][c]->SetStats(kFALSE);
2627    fEbECumulantsPro[cs][c]->SetMarkerColor(markerColor[cs]);
2628    fEbECumulantsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
2629    fEbECumulantsList->Add(fEbECumulantsPro[cs][c]);
2630   } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
2631  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2632  // Set all bin labels: TBI this can be implemented better, most likely...
2633  Int_t binNo[8]; for(Int_t c=0;c<fMaxCorrelator;c++){binNo[c]=1;} // TBI hardwired 8, shall be fMaxCorrelator
2634  for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++) 
2635  {
2636   if(fSkipZeroHarmonics && 0==n1){continue;}
2637   if(fCalculateAll)
2638   {
2639    fEbECumulantsPro[0][0]->GetXaxis()->SetBinLabel(binNo[0],Form("Cos(%d)",n1));
2640    fEbECumulantsPro[1][0]->GetXaxis()->SetBinLabel(binNo[0]++,Form("Sin(%d)",n1));
2641    nToBeFilled[0]++; 
2642   }
2643   if(1==fDontGoBeyond){continue;}
2644   for(Int_t n2=n1;n2<=fMaxHarmonic;n2++) 
2645   {
2646    if(fSkipZeroHarmonics && 0==n2){continue;}
2647    if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2)) 
2648       || (fCalculateSameIsotropic && 0==n1+n2 &&  TMath::Abs(n1)==TMath::Abs(n2)))
2649    {  
2650     fEbECumulantsPro[0][1]->GetXaxis()->SetBinLabel(binNo[1],Form("Cos(%d,%d)",n1,n2));
2651     fEbECumulantsPro[1][1]->GetXaxis()->SetBinLabel(binNo[1]++,Form("Sin(%d,%d)",n1,n2));
2652     nToBeFilled[1]++; 
2653    }
2654    if(2==fDontGoBeyond){continue;}
2655    for(Int_t n3=n2;n3<=fMaxHarmonic;n3++) 
2656    {
2657     if(fSkipZeroHarmonics && 0==n3){continue;}
2658     if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3))
2659        || (fCalculateSameIsotropic && 0==n1+n2+n3 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)))
2660     {  
2661      fEbECumulantsPro[0][2]->GetXaxis()->SetBinLabel(binNo[2],Form("Cos(%d,%d,%d)",n1,n2,n3));
2662      fEbECumulantsPro[1][2]->GetXaxis()->SetBinLabel(binNo[2]++,Form("Sin(%d,%d,%d)",n1,n2,n3));
2663      nToBeFilled[2]++; 
2664     }
2665     if(3==fDontGoBeyond){continue;}
2666     for(Int_t n4=n3;n4<=fMaxHarmonic;n4++) 
2667     {
2668      if(fSkipZeroHarmonics && 0==n4){continue;}
2669      if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4))
2670        || (fCalculateSameIsotropic && 0==n1+n2+n3+n4 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)))
2671      {   
2672       fEbECumulantsPro[0][3]->GetXaxis()->SetBinLabel(binNo[3],Form("Cos(%d,%d,%d,%d)",n1,n2,n3,n4));
2673       fEbECumulantsPro[1][3]->GetXaxis()->SetBinLabel(binNo[3]++,Form("Sin(%d,%d,%d,%d)",n1,n2,n3,n4)); 
2674       nToBeFilled[3]++; 
2675      } 
2676      if(4==fDontGoBeyond){continue;}
2677      for(Int_t n5=n4;n5<=fMaxHarmonic;n5++) 
2678      {
2679       if(fSkipZeroHarmonics && 0==n5){continue;}
2680       if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5) 
2681          || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5))
2682          || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5)))
2683       {   
2684        fEbECumulantsPro[0][4]->GetXaxis()->SetBinLabel(binNo[4],Form("Cos(%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5));
2685        fEbECumulantsPro[1][4]->GetXaxis()->SetBinLabel(binNo[4]++,Form("Sin(%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5));
2686        nToBeFilled[4]++; 
2687       }
2688       if(5==fDontGoBeyond){continue;}
2689       for(Int_t n6=n5;n6<=fMaxHarmonic;n6++) 
2690       {
2691        if(fSkipZeroHarmonics && 0==n6){continue;}
2692        if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6)  
2693           || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2694               && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6))
2695           || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
2696               && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6)))
2697        {   
2698         fEbECumulantsPro[0][5]->GetXaxis()->SetBinLabel(binNo[5],Form("Cos(%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6));
2699         fEbECumulantsPro[1][5]->GetXaxis()->SetBinLabel(binNo[5]++,Form("Sin(%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6));
2700         nToBeFilled[5]++; 
2701        }
2702        if(6==fDontGoBeyond){continue;}
2703        for(Int_t n7=n6;n7<=fMaxHarmonic;n7++) 
2704        {
2705         if(fSkipZeroHarmonics && 0==n7){continue;}
2706         if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7) 
2707            || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2708                && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7))
2709            || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2710                && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7)))
2711         {   
2712          fEbECumulantsPro[0][6]->GetXaxis()->SetBinLabel(binNo[6],Form("Cos(%d,%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6,n7));
2713          fEbECumulantsPro[1][6]->GetXaxis()->SetBinLabel(binNo[6]++,Form("Sin(%d,%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6,n7));
2714          nToBeFilled[6]++; 
2715         }
2716         if(7==fDontGoBeyond){continue;}
2717         for(Int_t n8=n7;n8<=fMaxHarmonic;n8++) 
2718         {
2719          if(fSkipZeroHarmonics && 0==n8){continue;}
2720          if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8) 
2721             || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) 
2722                 && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7) && TMath::Abs(n1)==TMath::Abs(n8))
2723             || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) 
2724                 && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7) 
2725                 && TMath::Abs(n1)==TMath::Abs(n8)))
2726          {    
2727           fEbECumulantsPro[0][7]->GetXaxis()->SetBinLabel(binNo[7],Form("Cos(%d,%d,%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6,n7,n8));
2728           fEbECumulantsPro[1][7]->GetXaxis()->SetBinLabel(binNo[7]++,Form("Sin(%d,%d,%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6,n7,n8));
2729           nToBeFilled[7]++; 
2730          }
2731         } // for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
2732        } // for(Int_t n7=n6;n7<=fMaxHarmonic;n7++) 
2733       } // for(Int_t n6=n5;n6<=fMaxHarmonic;n6++) 
2734      } // for(Int_t n5=n4;n5<=fMaxHarmonic;n5++) 
2735     } // for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)   
2736    } // for(Int_t n3=n2;n3<=fMaxHarmonic;n3++) 
2737   } // for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
2738  } // for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++) 
2739
2740  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2741  {
2742   for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2743   {
2744    fEbECumulantsPro[cs][c]->SetTitle(Form("%d-p e-b-e cumulants, %s terms, %d/%d in total",c+1,sCosSin[cs].Data(),nToBeFilled[c],nBinsTitle[c]));
2745    fEbECumulantsPro[cs][c]->GetXaxis()->SetRangeUser(0.,fEbECumulantsPro[cs][c]->GetBinLowEdge(nToBeFilled[c]+1));
2746   }
2747  } 
2748  
2749 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForEbECumulants()
2750
2751 //=======================================================================================================================
2752
2753 void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForNestedLoops()
2754 {
2755  // Book all the stuff for nested loops.
2756
2757  // TBI this method is just ugly, who implemented it like this... 
2758
2759  // a) Set default harmonic values; 
2760  // b) Book the profile holding all the flags for nested loops;
2761  // c) Book the profile holding all results for nested loops (cosine);
2762  // d) Book the profile holding all results for nested loops (sinus);
2763  // e) Book the profile holding all results for differential nested loops.
2764
2765  // a) Set default harmonic values: 
2766  //delete gRandom; // TBI this is not safe here, 
2767  //gRandom = new TRandom3(0);
2768  Int_t h1 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1); // TBI reimplement all these lines eventually 
2769  Int_t h2 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
2770  Int_t h3 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
2771  Int_t h4 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
2772  Int_t h5 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
2773  Int_t h6 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
2774  Int_t h7 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
2775  Int_t h8 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
2776
2777  // REMARK: This values can be overriden in a steering macro via 
2778  // mpc->GetNestedLoopsFlagsPro()->SetBinContent(<binNo>,<value>);
2779
2780  // b) Book the profile holding all the flags for nested loops:
2781  fNestedLoopsFlagsPro = new TProfile("fNestedLoopsFlagsPro","Flags for nested loops",10,0,10);
2782  fNestedLoopsFlagsPro->SetTickLength(-0.01,"Y");
2783  fNestedLoopsFlagsPro->SetMarkerStyle(25);
2784  fNestedLoopsFlagsPro->SetLabelSize(0.03);
2785  fNestedLoopsFlagsPro->SetLabelOffset(0.02,"Y");
2786  fNestedLoopsFlagsPro->SetStats(kFALSE);
2787  fNestedLoopsFlagsPro->SetFillColor(kGray);
2788  fNestedLoopsFlagsPro->SetLineColor(kBlack);
2789  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(1,"fCrossCheckWithNestedLoops"); fNestedLoopsFlagsPro->Fill(0.5,fCrossCheckWithNestedLoops);
2790  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(2,"h_{1}"); fNestedLoopsFlagsPro->Fill(1.5,h1);
2791  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(3,"h_{2}"); fNestedLoopsFlagsPro->Fill(2.5,h2);
2792  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(4,"h_{3}"); fNestedLoopsFlagsPro->Fill(3.5,h3);
2793  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(5,"h_{4}"); fNestedLoopsFlagsPro->Fill(4.5,h4);
2794  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(6,"h_{5}"); fNestedLoopsFlagsPro->Fill(5.5,h5);
2795  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(7,"h_{6}"); fNestedLoopsFlagsPro->Fill(6.5,h6);
2796  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(8,"h_{7}"); fNestedLoopsFlagsPro->Fill(7.5,h7);
2797  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(9,"h_{8}"); fNestedLoopsFlagsPro->Fill(8.5,h8);
2798  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(10,"fCrossCheckDiffWithNestedLoops"); fNestedLoopsFlagsPro->Fill(9.5,fCrossCheckDiffWithNestedLoops);
2799  fNestedLoopsList->Add(fNestedLoopsFlagsPro);
2800
2801  // c) Book the profile holding all results for nested loops (cosine):
2802  fNestedLoopsResultsCosPro = new TProfile("fNestedLoopsResultsCosPro","Nested loops results (cosine)",16,0,16);
2803  fNestedLoopsResultsCosPro->SetTickLength(-0.01,"Y");
2804  fNestedLoopsResultsCosPro->SetMarkerStyle(25);
2805  fNestedLoopsResultsCosPro->SetLabelSize(0.02);
2806  fNestedLoopsResultsCosPro->SetLabelOffset(0.02,"Y");
2807  fNestedLoopsResultsCosPro->SetStats(kFALSE);
2808  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(1,Form("N: 1p(%d)",h1));
2809  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(2,Form("Q: 1p(%d)",h1));
2810  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(3,Form("N: 2p(%d,%d)",h1,h2));
2811  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(4,Form("Q: 2p(%d,%d)",h1,h2));
2812  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(5,Form("N: 3p(%d,%d,%d)",h1,h2,h3));
2813  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(6,Form("Q: 3p(%d,%d,%d)",h1,h2,h3));
2814  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(7,Form("N: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
2815  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(8,Form("Q: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
2816  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(9,Form("N: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
2817  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(10,Form("Q: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
2818  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(11,Form("N: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
2819  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(12,Form("Q: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
2820  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(13,Form("N: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
2821  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(14,Form("Q: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
2822  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(15,Form("N: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
2823  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(16,Form("Q: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
2824  if(fCrossCheckWithNestedLoops){fNestedLoopsList->Add(fNestedLoopsResultsCosPro);} else{delete fNestedLoopsResultsCosPro;}
2825
2826  // d) Book the profile holding all results for nested loops (sinus):
2827  fNestedLoopsResultsSinPro = new TProfile("fNestedLoopsResultsSinPro","Nested loops results (sinus)",16,0,16);
2828  fNestedLoopsResultsSinPro->SetTickLength(-0.01,"Y");
2829  fNestedLoopsResultsSinPro->SetMarkerStyle(25);
2830  fNestedLoopsResultsSinPro->SetLabelSize(0.02);
2831  fNestedLoopsResultsSinPro->SetLabelOffset(0.02,"Y");
2832  fNestedLoopsResultsSinPro->SetStats(kFALSE);
2833  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(1,Form("N: 1p(%d)",h1));
2834  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(2,Form("Q: 1p(%d)",h1));
2835  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(3,Form("N: 2p(%d,%d)",h1,h2));
2836  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(4,Form("Q: 2p(%d,%d)",h1,h2));
2837  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(5,Form("N: 3p(%d,%d,%d)",h1,h2,h3));
2838  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(6,Form("Q: 3p(%d,%d,%d)",h1,h2,h3));
2839  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(7,Form("N: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
2840  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(8,Form("Q: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
2841  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(9,Form("N: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
2842  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(10,Form("Q: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
2843  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(11,Form("N: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
2844  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(12,Form("Q: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
2845  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(13,Form("N: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
2846  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(14,Form("Q: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
2847  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(15,Form("N: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
2848  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(16,Form("Q: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
2849  if(fCrossCheckWithNestedLoops){fNestedLoopsList->Add(fNestedLoopsResultsSinPro);} else{delete fNestedLoopsResultsSinPro;}
2850
2851  // e) Book the profile holding all results for differential nested loops:
2852  fNestedLoopsDiffResultsPro = new TProfile("fNestedLoopsDiffResultsPro","Differential nested loops results",1,0.,1.);
2853  fNestedLoopsDiffResultsPro->SetStats(kFALSE);
2854  if(fCrossCheckDiffWithNestedLoops){fNestedLoopsList->Add(fNestedLoopsDiffResultsPro);} else{delete fNestedLoopsDiffResultsPro;}
2855
2856 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForNestedLoops()
2857
2858 //=======================================================================================================================
2859
2860 void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForStandardCandles()
2861 {
2862  // Book all the stuff for 'standard candles'.
2863
2864  // a) Book the profile holding all the flags for 'standard candles';
2865  // b) Book the histogram holding all results for 'standard candles';
2866  // c) Book TProfile2D *fProductsSCPro. 
2867
2868  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForStandardCandles()";
2869
2870  // a) Book the profile holding all the flags for 'standard candles':
2871  fStandardCandlesFlagsPro = new TProfile("fStandardCandlesFlagsPro","Flags for standard candles",2,0,2);
2872  fStandardCandlesFlagsPro->SetTickLength(-0.01,"Y");
2873  fStandardCandlesFlagsPro->SetMarkerStyle(25);
2874  fStandardCandlesFlagsPro->SetLabelSize(0.03);
2875  fStandardCandlesFlagsPro->SetLabelOffset(0.02,"Y");
2876  fStandardCandlesFlagsPro->SetStats(kFALSE);
2877  fStandardCandlesFlagsPro->SetFillColor(kGray);
2878  fStandardCandlesFlagsPro->SetLineColor(kBlack);
2879  fStandardCandlesFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateStandardCandles"); fStandardCandlesFlagsPro->Fill(0.5,fCalculateStandardCandles);
2880  fStandardCandlesFlagsPro->GetXaxis()->SetBinLabel(2,"fPropagateErrorSC"); fStandardCandlesFlagsPro->Fill(1.5,fPropagateErrorSC);
2881  fStandardCandlesList->Add(fStandardCandlesFlagsPro);
2882
2883  if(!fCalculateStandardCandles){return;} // TBI is this safe like this? 
2884
2885  // b) Book the histogram holding all results for 'standard candles':
2886  Int_t nBins = fMaxHarmonic*(fMaxHarmonic-1)/2;
2887  fStandardCandlesHist = new TH1D("fStandardCandlesHist","Standard candles (SC)",nBins,0.,1.*nBins); 
2888  fStandardCandlesHist->SetStats(kFALSE);
2889  fStandardCandlesHist->SetMarkerStyle(kFullSquare);
2890  fStandardCandlesHist->SetMarkerColor(kBlue);
2891  fStandardCandlesHist->SetLineColor(kBlue);
2892  Int_t binNo = 1;
2893  for(Int_t n1=-fMaxHarmonic;n1<=-2;n1++)
2894  {
2895   for(Int_t n2=n1+1;n2<=-1;n2++)
2896   {
2897    fStandardCandlesHist->GetXaxis()->SetBinLabel(binNo++,Form("SC(%d,%d,%d,%d)",n1,n2,-1*n2,-1*n1));
2898   }
2899  }
2900  if(binNo-1 != nBins){Fatal(sMethodName.Data(),"Well, binNo-1 != nBins ... :'( ");}
2901  fStandardCandlesList->Add(fStandardCandlesHist);
2902
2903  if(!fPropagateErrorSC){return;} // TBI is this safe like this? 
2904
2905  // c) Book TProfile2D *fProductsSCPro:
2906  Int_t nBins2D = fMaxHarmonic + fMaxHarmonic*(fMaxHarmonic-1)/2; // #2-p + #4-p distinct correlators in SC context 
2907  if(nBins2D<=0){Fatal(sMethodName.Data(),"nBins2D<=0");} // well, just in case...
2908  fProductsSCPro = new TProfile2D("fProductsSCPro","Products of correlations",nBins2D,0.,nBins2D,nBins2D,0.,nBins2D);
2909  fProductsSCPro->SetStats(kFALSE);
2910  fProductsSCPro->Sumw2();
2911  for(Int_t b=1;b<=fMaxHarmonic;b++) // 2-p correlators
2912  {
2913   fProductsSCPro->GetXaxis()->SetBinLabel(b,Form("Cos(%d,%d)",-(fMaxHarmonic+1-b),(fMaxHarmonic+1-b)));
2914   fProductsSCPro->GetYaxis()->SetBinLabel(b,Form("Cos(%d,%d)",-(fMaxHarmonic+1-b),(fMaxHarmonic+1-b)));
2915  } // for(Int_t b=1;b<=fMaxHarmonic;b++) // 2-p correlators
2916  for(Int_t b=fMaxHarmonic+1;b<=nBins2D;b++) // 4-p correlators
2917  {
2918   TString sBinLabel = fStandardCandlesHist->GetXaxis()->GetBinLabel(b-fMaxHarmonic);
2919   if(sBinLabel.EqualTo("")){Fatal(sMethodName.Data(),"sBinLabel.EqualTo...");}
2920   fProductsSCPro->GetXaxis()->SetBinLabel(b,sBinLabel.ReplaceAll("SC","Cos").Data());
2921   fProductsSCPro->GetYaxis()->SetBinLabel(b,sBinLabel.ReplaceAll("SC","Cos").Data());
2922  } // for(Int_t b=fMaxHarmonic+1;b<=nBins2D;b++) // 4-p correlators
2923  fStandardCandlesList->Add(fProductsSCPro);
2924
2925 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForStandardCandles()
2926
2927 //=======================================================================================================================
2928
2929 void AliFlowAnalysisWithMultiparticleCorrelations::GetOutputHistograms(TList *histList)
2930 {
2931  // Get pointers for everything and everywhere from the base list "fHistList". 
2932
2933  // a) Get pointer for base list fHistList;
2934  // b) Get pointer for profile holding internal flags and, well, set again all flags;
2935  // c) Get pointers for control histograms;
2936  // d) Get pointers for Q-vector;
2937  // e) Get pointers for correlations;
2938  // f) Get pointers for 'standard candles';
2939  // g) Get pointers for Q-cumulants;
2940  // h) Get pointers for differential correlations.
2941
2942  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetOutputHistograms(TList *histList)";
2943
2944  // a) Get pointer for base list fHistList and profile holding internal flags;
2945  fHistList = histList; 
2946  if(!fHistList){Fatal(sMethodName.Data(),"fHistList is malicious today...");}
2947
2948  // b) Get pointer for profile holding internal flags and, well, set again all flags:
2949  fInternalFlagsPro = dynamic_cast<TProfile*>(fHistList->FindObject("fInternalFlagsPro"));
2950  if(!fInternalFlagsPro){Fatal(sMethodName.Data(),"fInternalFlagsPro");}
2951  fUseInternalFlags = fInternalFlagsPro->GetBinContent(1);
2952  fMinNoRPs = (Int_t)fInternalFlagsPro->GetBinContent(2);
2953  fMaxNoRPs = (Int_t)fInternalFlagsPro->GetBinContent(3);
2954  fExactNoRPs = (Int_t)fInternalFlagsPro->GetBinContent(4);
2955  fPropagateError = (Bool_t)fInternalFlagsPro->GetBinContent(5);
2956
2957  // c) Get pointers for control histograms:
2958  this->GetPointersForControlHistograms(); 
2959
2960  // d) Get pointers for Q-vector:
2961  this->GetPointersForQvector(); 
2962
2963  // e) Get pointers for correlations:
2964  this->GetPointersForCorrelations(); 
2965
2966  // f) Get pointers for 'standard candles':
2967  this->GetPointersForStandardCandles(); 
2968
2969  // g) Get pointers for Q-cumulants:
2970  this->GetPointersForQcumulants(); 
2971
2972  // h) Get pointers for differential correlations:
2973  this->GetPointersForDiffCorrelations(); 
2974   
2975 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetOutputHistograms(TList *histList)
2976
2977 //=======================================================================================================================
2978
2979 void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForQvector()
2980 {
2981  // Get pointers for Q-vector objects.
2982
2983  // a) Get pointer for fQvectorList; 
2984  // b) Get pointer for fQvectorFlagsPro;
2985  // c) Set again all flags.
2986
2987  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForQvector()";
2988
2989  // a) Get pointer for fQvectorList: 
2990  fQvectorList = dynamic_cast<TList*>(fHistList->FindObject("Q-vectors")); 
2991  if(!fQvectorList){Fatal(sMethodName.Data(),"fQvectorList");}
2992
2993  // b) Get pointer for fQvectorFlagsPro: 
2994  fQvectorFlagsPro = dynamic_cast<TProfile*>(fQvectorList->FindObject("fQvectorFlagsPro"));
2995  if(!fQvectorFlagsPro){Fatal(sMethodName.Data(),"fQvectorFlagsPro");}
2996
2997  // c) Set again all flags:
2998  fCalculateQvector = (Bool_t)fQvectorFlagsPro->GetBinContent(1);
2999  fCalculateDiffQvectors = (Bool_t)fQvectorFlagsPro->GetBinContent(2);
3000
3001 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForQvector()
3002
3003 //=======================================================================================================================
3004
3005 void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForStandardCandles()
3006 {
3007  // Get pointers for 'standard candles'.
3008
3009  // a) Get pointer for fStandardCandlesList; 
3010  // b) Get pointer for fStandardCandlesFlagsPro;
3011  // c) Set again all flags; 
3012  // d) Get pointer TH1D *fStandardCandlesHist;
3013  // e) Get pointer TProfile2D *fProductsSCPro.
3014
3015  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForStandardCandles()";
3016
3017  // a) Get pointer for fStandardCandlesList: 
3018  fStandardCandlesList = dynamic_cast<TList*>(fHistList->FindObject("Standard Candles"));
3019  if(!fStandardCandlesList){Fatal(sMethodName.Data(),"fStandardCandlesList");}
3020
3021  // b) Get pointer for fStandardCandlesFlagsPro: 
3022  fStandardCandlesFlagsPro = dynamic_cast<TProfile*>(fStandardCandlesList->FindObject("fStandardCandlesFlagsPro"));
3023  if(!fStandardCandlesFlagsPro){Fatal(sMethodName.Data(),"fStandardCandlesFlagsPro");}
3024
3025  // c) Set again all flags: 
3026  fCalculateStandardCandles = (Bool_t)fStandardCandlesFlagsPro->GetBinContent(1);
3027  fPropagateErrorSC = (Bool_t)fStandardCandlesFlagsPro->GetBinContent(2);
3028
3029  if(!fCalculateStandardCandles){return;} // TBI is this safe enough
3030
3031  // d) Get pointer TH1D *fStandardCandlesHist: 
3032  fStandardCandlesHist = dynamic_cast<TH1D*>(fStandardCandlesList->FindObject("fStandardCandlesHist"));
3033
3034  if(!fPropagateErrorSC){return;} // TBI is this safe enough
3035
3036  // e) Get pointer TProfile2D *fProductsSCPro:
3037  fProductsSCPro = dynamic_cast<TProfile2D*>(fStandardCandlesList->FindObject("fProductsSCPro"));
3038  
3039 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForStandardCandles()
3040
3041 //=======================================================================================================================
3042
3043 void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForControlHistograms()
3044 {
3045  // Get pointers for control histograms.
3046
3047  // a) Get pointer for fControlHistogramsList; TBI
3048  // b) Get pointer for fControlHistogramsFlagsPro; TBI
3049  // c) Set again all flags; TBI
3050  // d) Get pointers to TH1D *fKinematicsHist[2][3]; TBI
3051  // e) Get pointers to TH1D *fMultDistributionsHist[3]; TBI
3052  // f) Get pointers to TH2D *fMultCorrelationsHist[3]. TBI
3053
3054  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForControlHistograms()";
3055
3056  // a) Get pointer for fControlHistogramsList: TBI
3057  fControlHistogramsList = dynamic_cast<TList*>(fHistList->FindObject("Control Histograms"));
3058  if(!fControlHistogramsList){Fatal(sMethodName.Data(),"fControlHistogramsList");}
3059
3060  // b) Get pointer for fControlHistogramsFlagsPro: TBI
3061  fControlHistogramsFlagsPro = dynamic_cast<TProfile*>(fControlHistogramsList->FindObject("fControlHistogramsFlagsPro"));
3062  if(!fControlHistogramsFlagsPro){Fatal(sMethodName.Data(),"fControlHistogramsFlagsPro");}
3063
3064  // c) Set again all flags:
3065  fFillControlHistograms = fControlHistogramsFlagsPro->GetBinContent(1);
3066  fFillKinematicsHist = fControlHistogramsFlagsPro->GetBinContent(2);
3067  fFillMultDistributionsHist = fControlHistogramsFlagsPro->GetBinContent(3);
3068  fFillMultCorrelationsHist = fControlHistogramsFlagsPro->GetBinContent(4);
3069
3070  if(!fFillControlHistograms){return;} // TBI