AliPhysics  608b256 (608b256)
AliFlowAnalysisWithMultiparticleCorrelations.cxx
Go to the documentation of this file.
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 
27 
28 using std::endl;
29 using std::cout;
30 using std::flush;
31 using std::ofstream;
32 using std::ios;
33 
34 //================================================================================================================
35 
37 
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  fSelectRandomlyRPs(kFALSE),
51  fnSelectedRandomlyRPs(-44),
52  fRandomIndicesRPs(NULL),
53  // 1.) Control histograms:
54  fControlHistogramsList(NULL),
55  fControlHistogramsFlagsPro(NULL),
56  fFillControlHistograms(kFALSE),
57  fFillKinematicsHist(kFALSE),
58  fFillMultDistributionsHist(kFALSE),
59  fFillMultCorrelationsHist(kFALSE),
60  fSkipSomeIntervals(kFALSE),
61  fNumberOfSkippedRPParticles(0),
62  // 2.) Q-vector:
63  fQvectorList(NULL),
64  fQvectorFlagsPro(NULL),
65  fCalculateQvector(kFALSE),
66  fCalculateDiffQvectors(kFALSE),
67  // 3.) Correlations:
68  fCorrelationsList(NULL),
69  fCorrelationsFlagsPro(NULL),
70  fCalculateCorrelations(kFALSE),
71  fMaxHarmonic(6), // TBI this shall not be hardwired in the ideal world...
72  fMaxCorrelator(8),
73  fCalculateIsotropic(kFALSE),
74  fCalculateSame(kFALSE),
75  fSkipZeroHarmonics(kFALSE),
76  fCalculateSameIsotropic(kFALSE),
77  fCalculateAll(kFALSE),
78  fDontGoBeyond(0),
79  fCalculateOnlyForHarmonicQC(kFALSE),
80  fCalculateOnlyForSC(kFALSE),
81  fCalculateOnlyCos(kFALSE),
82  fCalculateOnlySin(kFALSE),
83  // 4.) Event-by-event cumulants:
84  fEbECumulantsList(NULL),
85  fEbECumulantsFlagsPro(NULL),
86  fCalculateEbECumulants(kFALSE),
87  // 5.) Weights:
88  fWeightsList(NULL),
89  fWeightsFlagsPro(NULL),
90  // 6.) Nested loops:
91  fNestedLoopsList(NULL),
92  fNestedLoopsFlagsPro(NULL),
93  fCrossCheckWithNestedLoops(kFALSE),
94  fCrossCheckDiffWithNestedLoops(kFALSE),
95  fNestedLoopsResultsCosPro(NULL),
96  fNestedLoopsResultsSinPro(NULL),
97  fNestedLoopsDiffResultsPro(NULL),
98  // 7.) 'Standard candles':
99  fStandardCandlesList(NULL),
100  fStandardCandlesFlagsPro(NULL),
101  fCalculateStandardCandles(kFALSE),
102  fPropagateErrorSC(kTRUE),
103  fStandardCandlesHist(NULL),
104  fProductsSCPro(NULL),
105  // 8.) Q-cumulants:
106  fQcumulantsList(NULL),
107  fQcumulantsFlagsPro(NULL),
108  fCalculateQcumulants(kFALSE),
109  fHarmonicQC(2),
110  fPropagateErrorQC(kTRUE),
111  fQcumulantsHist(NULL),
112  fReferenceFlowHist(NULL),
113  fProductsQCPro(NULL),
114  // 9.) Differential correlations:
115  fDiffCorrelationsList(NULL),
116  fDiffCorrelationsFlagsPro(NULL),
117  fCalculateDiffCorrelations(kFALSE),
118  fCalculateDiffCos(kTRUE),
119  fCalculateDiffSin(kFALSE),
120  fCalculateDiffCorrelationsVsPt(kTRUE),
121  fUseDefaultBinning(kTRUE),
122  fnDiffBins(-44),
123  fRangesDiffBins(NULL),
124  fDiffBinNo(-1),
125  // 10.) Symmetry plane correlations:
126  fSymmetryPlanesList(NULL),
127  fSymmetryPlanesFlagsPro(NULL),
128  fCalculateSymmetryPlanes(kFALSE),
129  // 11.) Eta gaps:
130  fEtaGapsList(NULL),
131  fEtaGapsFlagsPro(NULL),
132  fCalculateEtaGaps(kFALSE),
133  fLowestHarmonicEtaGaps(1),
134  fHighestHarmonicEtaGaps(6)
135  {
136  // Constructor.
137 
138  // a) Book grandmother of all lists;
139  // b) Initialize all arrays.
140 
141  // a) Book grandmother of all lists:
142  fHistList = new TList();
143  fHistList->SetName("cobjMPC");
144  fHistList->SetOwner(kTRUE);
145 
146  // b) Initialize all arrays:
147  this->InitializeArraysForControlHistograms();
148  this->InitializeArraysForQvector();
149  this->InitializeArraysForCorrelations();
150  this->InitializeArraysForEbECumulants();
151  this->InitializeArraysForWeights();
152  this->InitializeArraysForQcumulants();
153  this->InitializeArraysForDiffCorrelations();
154  this->InitializeArraysForSymmetryPlanes();
155  this->InitializeArraysForNestedLoops();
156  this->InitializeArraysForEtaGaps();
157  // TBI test for the workflow
158 
159  // c) Determine seed for gRandom:
160  delete gRandom;
161  gRandom = new TRandom3(0); // since 0 is in the argument, the seed is determined uniquely in space and time via TUUID
162  // TBI synchronize this eventually with seed set 'on-the-fly'
163 
164  } // end of AliFlowAnalysisWithMultiparticleCorrelations::AliFlowAnalysisWithMultiparticleCorrelations()
165 
166 //================================================================================================================
167 
169 {
170  // Destructor.
171 
172  delete fHistList;
173 
174 } // end of AliFlowAnalysisWithMultiparticleCorrelations::~AliFlowAnalysisWithMultiparticleCorrelations()
175 
176 //================================================================================================================
177 
179 {
180  // Well, this is method Init().
181 
182  // a) Trick to avoid name clashes, part 1;
183  // b) Cross-check the initial settings before starting this adventure;
184  // c) Book all objects;
185  // d) Set all flags;
186  // *) Trick to avoid name clashes, part 2.
187 
188  // a) Trick to avoid name clashes, part 1:
189  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
190  TH1::AddDirectory(kFALSE);
191 
192  // b) Cross-check the initial settings before starting this adventure:
193  this->CrossCheckSettings();
194 
195  // c) Book all objects:
196  this->BookAndNestAllLists();
197  this->BookEverythingForBase();
199  this->BookEverythingForQvector();
200  this->BookEverythingForWeights();
208  this->BookEverythingForEtaGaps();
209 
210  // d) Set all flags:
211  // ...
212 
213  // *) Trick to avoid name clashes, part 2:
214  TH1::AddDirectory(oldHistAddStatus);
215 
216 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::Init()
217 
218 //================================================================================================================
219 
221 {
222  // Running over data only in this method.
223 
224  // a) Cross-check internal flags;
225  // b) Cross-check all pointers used in this method;
226  // c) Determine random indices;
227  // d) Fill control histograms;
228  // e) Fill Q-vector components;
229  // f) Calculate multi-particle correlations from Q-vector components;
230  // g) Calculate e-b-e cumulants;
231  // h) Calculate symmetry plane correlations;
232  // i) Calculate 2-p correlations with eta gaps;
233  // j) Reset Q-vector components;
234  // k) Cross-check results with nested loops;
235  // l) Dump the points;
236  // m) Reset array holding shuffled indices for RPs.
237 
238  // a) Cross-check internal flags:
239  if(fUseInternalFlags){if(!this->CrossCheckInternalFlags(anEvent)){return;}}
240 
241  // b) Cross-check all pointers used in this method:
242  this->CrossCheckPointersUsedInMake(); // TBI shall I call this method first
243 
244  // c) Determine random indices:
246  {
247  if(anEvent->GetNumberOfRPs() < fnSelectedRandomlyRPs){return;}
248  this->DetermineRandomIndices(anEvent);
249  if(!fRandomIndicesRPs){return;}
250  } // TBI hw RPs in flag, make it more general
251 
252  // TBI temp gym: Remove this code eventually
254  {
256  Int_t nTracks = anEvent->NumberOfTracks(); // TBI shall I promote this to data member?
257  for(Int_t t=0;t<nTracks;t++) // loop over all tracks
258  {
259  AliFlowTrackSimple *pTrack = anEvent->GetTrack(t);
260  if(!pTrack){printf("\n AAAARGH: pTrack is NULL in MPC::FCH() !!!!");continue;}
261  if(!pTrack->InRPSelection()){continue;}
262  if(pTrack)
263  {
264  Double_t dPhi = pTrack->Phi();
265  //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
266  //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
267  Double_t dPt = pTrack->Pt();
268  Double_t dEta = pTrack->Eta();
269  Double_t dPhiPtEta[3] = {dPhi,dPt,dEta};
270 
271  // Skip some intervals: TBI promote eventually to AFTC class
272  Bool_t bPasses = kTRUE;
273  Bool_t bAlreadyCounted = kFALSE;
274  for(Int_t ppe=0;ppe<3;ppe++)
275  {
276  if(!bPasses){break;} // found one kinematic window which shall be skipped
277  for(Int_t b=0;b<10;b+=2)
278  {
279  if(-44==(Int_t)fSkip[ppe][b]){continue;}
280  if(dPhiPtEta[ppe]>=fSkip[ppe][b] && dPhiPtEta[ppe]<fSkip[ppe][b+1])
281  {
282  bPasses = kFALSE;
283  if(!bAlreadyCounted)
284  {
286  bAlreadyCounted = kTRUE;
287  } // if(bAlreadyCounted)
288  break;
289  } // TBI this is a clear bug when this setter is used multiple times...
290  } // for(Int_t b=0;b<10;b++)
291  } // for(Int_t ppe=0;ppe<3;ppe++)
292  if(!bPasses){continue;}
293 
294  } // if(pTrack)
295  } // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
296  } // if(fSkipSomeIntervals)
297  if(fSkipSomeIntervals) // TBI tmp gym
298  {
299  if(anEvent->GetNumberOfRPs() - fNumberOfSkippedRPParticles < 8){return;} // TBI tmp gym
300  }
301 
302 
303 
304  // d) Fill control histograms:
305  if(fFillControlHistograms){this->FillControlHistograms(anEvent);}
306 
307  // e) Fill Q-vector components:
309 
310  // f) Calculate multi-particle correlations from Q-vector components:
311  if(fCalculateCorrelations){this->CalculateCorrelations(anEvent);}
313 
314  // g) Calculate e-b-e cumulants:
315  if(fCalculateEbECumulants){this->CalculateEbECumulants(anEvent);}
316 
317  // h) Calculate symmetry plane correlations:
319 
320  // i) Calculate 2-p correlations with eta gaps:
321  if(fCalculateEtaGaps){this->CalculateEtaGaps(anEvent);}
322 
323  // j) Reset Q-vector components:
325 
326  // k) Cross-check results with nested loops:
329 
330  // l) Dump the points:
331  if(fDumpThePoints){this->DumpThePoints(anEvent);}
332 
333  // m) Reset array holding shuffled indices for RPs:
335 
336 } // end of AliFlowAnalysisWithMultiparticleCorrelations::Make(AliFlowEventSimple *anEvent)
337 
338 //=======================================================================================================================
339 
341 {
342  // Closing the curtains.
343 
344  // a) Cross-check pointers used in this method;
345  // b) Calculate 'standard candles';
346  // c) Calculate Q-cumulants.
347 
348  // a) Cross-check pointers used in this method:
350 
351  // b) Calculate 'standard candles':
353 
354  // c) Calculate Q-cumulants:
356 
357  // ...
358 
359  printf("\n ... Closing the curtains ... \n\n");
360 
361 } // end of AliFlowAnalysisWithMultiparticleCorrelations::Finish()
362 
363 //=======================================================================================================================
364 
366 {
367  // Cross-check all pointers used in method Finish().
368 
369  // a) Correlations;
370  // b) 'Standard candles';
371  // c) Q-cumulants.
372 
373  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInFinish()";
374 
375  // a) Correlations:
377  {
378  for(Int_t cs=0;cs<2;cs++)
379  {
380  if(fCalculateOnlyCos && 1==cs){continue;}
381  else if(fCalculateOnlySin && 0==cs){continue;}
382  for(Int_t c=0;c<fMaxCorrelator;c++)
383  {
384  if(c==fDontGoBeyond){continue;}
385  if(fCalculateOnlyForHarmonicQC && c%2==0){continue;}
386  if(!fCorrelationsPro[cs][c]){Fatal(sMethodName.Data(),"fCorrelationsPro[%d][%d] is NULL, for one reason or another...",cs,c);}
387  }
388  }
390  {Fatal(sMethodName.Data(),"fCalculateQcumulants && fPropagateErrorQC && !fProductsQCPro");}
391  } // if(fCalculateCorrelations)
392 
393  // b) 'Standard candles':
395  {
396  if(!fStandardCandlesHist){Fatal(sMethodName.Data(),"fStandardCandlesHist is NULL, for one reason or another...");}
398  {
399  if(!fProductsSCPro){Fatal(sMethodName.Data(),"fProductsSCPro is NULL, for one reason or another...");}
400  }
401  } // if(fCalculateStandardCandles)
402 
403  // c) Q-cumulants:
405  {
406  if(!fQcumulantsHist){Fatal(sMethodName.Data(),"fQcumulantsHist is NULL, for one reason or another...");}
407  if(!fReferenceFlowHist){Fatal(sMethodName.Data(),"fReferenceFlowHist is NULL, for one reason or another...");}
409  {
410  if(!fProductsQCPro){Fatal(sMethodName.Data(),"fProductsQCPro is NULL, for one reason or another...");}
411  }
412  } // if(fCalculateQcumulants)
413 
414 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInFinish()
415 
416 //=======================================================================================================================
417 
419 {
420  // Cross-check all pointers used in method Make().
421 
422  // a) Correlations;
423  // b) Event-by-event cumulants;
424  // c) ...
425 
426  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInMake()";
427 
428  // a) Correlations:
430  {
431  for(Int_t cs=0;cs<2;cs++)
432  {
433  if(fCalculateOnlyCos && 1==cs){continue;}
434  else if(fCalculateOnlySin && 0==cs){continue;}
435  for(Int_t c=0;c<fMaxCorrelator;c++)
436  {
437  if(c==fDontGoBeyond){continue;}
438  if(fCalculateOnlyForHarmonicQC && c%2==0){continue;}
439  if(!fCorrelationsPro[cs][c]){Fatal(sMethodName.Data(),"fCorrelationsPro[%d][%d] is NULL, for one reason or another...",cs,c);}
440  }
441  }
443  {Fatal(sMethodName.Data(),"fCalculateQcumulants && fPropagateErrorQC && !fProductsQCPro");}
444  } // if(fCalculateCorrelations)
445 
446  // b) Event-by-event cumulants:
448  {
449  for(Int_t cs=0;cs<2;cs++)
450  {
451  for(Int_t c=0;c<fMaxCorrelator;c++)
452  {
453  if(!fEbECumulantsPro[cs][c]){Fatal(sMethodName.Data(),"fEbECumulantsPro[%d][%d] is NULL, for one reason or another...",cs,c);}
454  }
455  }
456  } // if(fCalculateEbECumulants)
457 
458 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckPointersUsedInMake()
459 
460 //=======================================================================================================================
461 
463 {
464  // Determine random indices.
465 
466  // Fisher-Yates algorithm:
467  Int_t nPrim = anEvent->NumberOfTracks();
468  if(nPrim > 0)
469  {
470  fRandomIndicesRPs = new TArrayI(nPrim);
471  }
472  else
473  {
474  return;
475  }
476 
477  for(Int_t i=0;i<nPrim;i++)
478  {
479  fRandomIndicesRPs->AddAt(i,i);
480  }
481  for(Int_t i=nPrim-1;i>=1;i--)
482  {
483  Int_t j = gRandom->Integer(i+1);
484  Int_t temp = fRandomIndicesRPs->GetAt(j);
485  fRandomIndicesRPs->AddAt(fRandomIndicesRPs->GetAt(i),j);
486  fRandomIndicesRPs->AddAt(temp,i);
487  } // end of for(Int_t i=nPrim-1;i>=1;i--)
488 
489 } // void AliFlowAnalysisWithMultiparticleCorrelations::DetermineRandomIndices(AliFlowEventSimple *anEvent)
490 
491 //=======================================================================================================================
492 
494 {
495  // Calculate 'standard candles'.
496 
497  // 'Standard candle' (SC) is defined in terms of average (all-event!) correlations as follows:
498  // SC(-n1,-n2,n2,n1) = <<Cos(-n1,-n2,n2,n1)>> - <<Cos(-n1,n1)>>*<<Cos(-n2,n2)>>, n1 > n2.
499 
500  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateStandardCandles()";
501 
502  Int_t nBins = fStandardCandlesHist->GetNbinsX();
503  Int_t nBins2p = fCorrelationsPro[0][1]->GetNbinsX();
504  Int_t nBins4p = fCorrelationsPro[0][3]->GetNbinsX();
505  for(Int_t b=1;b<=nBins;b++)
506  {
507  // Standard candle:
508  Double_t dSCn1n2n2n1 = 0.; // SC(-n1,-n2,n2,n1)
509  Double_t dSCn1n2n2n1Err = 0.; // SC(-n1,-n2,n2,n1) stat. error
510  fPropagateError = kTRUE;
511 
512  // Correlations:
513  Double_t dCosn1n2n2n1 = 0.; // <<Cos(-n1,-n2,n2,n1)>>
514  Double_t dCosn1n2n2n1Err = 0.; // <<Cos(-n1,-n2,n2,n1)>> stat. error
515  Double_t dCosn1n1 = 0.; // <<Cos(-n1,n1)>>
516  Double_t dCosn1n1Err = 0.; // <<Cos(-n1,n1)>> stat. error
517  Double_t dCosn2n2 = 0.; // <<Cos(-n2,n2)>>
518  Double_t dCosn2n2Err = 0.; // <<Cos(-n2,n2)>> stat. error
519 
520  // Labels:
521  TString labeln1n2n2n1 = TString(fStandardCandlesHist->GetXaxis()->GetBinLabel(b)).ReplaceAll("SC","Cos");
522  TString n1 = TString(fStandardCandlesHist->GetXaxis()->GetBinLabel(b))(4);
523  TString n2 = TString(fStandardCandlesHist->GetXaxis()->GetBinLabel(b))(7);
524  if(n1.EqualTo("-") || n1.EqualTo(",")){Fatal(sMethodName.Data(),"n1.EqualTo...");}
525  if(n2.EqualTo("-") || n2.EqualTo(",")){Fatal(sMethodName.Data(),"n2.EqualTo...");}
526  TString labeln1n1 = Form("Cos(-%s,%s)",n1.Data(),n1.Data());
527  TString labeln2n2 = Form("Cos(-%s,%s)",n2.Data(),n2.Data());
528  //cout<<labeln1n2n2n1.Data()<<endl;
529  //cout<<labeln1n1.Data()<<endl;
530  //cout<<labeln2n2.Data()<<endl;
531  //cout<<endl;
532 
533  // Access <<Cos(-n1,-n2,n2,n1)>>:
534  for(Int_t b4p=1;b4p<=nBins4p;b4p++)
535  {
536  if(labeln1n2n2n1.EqualTo(fCorrelationsPro[0][3]->GetXaxis()->GetBinLabel(b4p)))
537  {
538  //cout<<labeln1n2n2n1.Data()<<endl;
539  dCosn1n2n2n1 = fCorrelationsPro[0][3]->GetBinContent(b4p);
540  dCosn1n2n2n1Err = fCorrelationsPro[0][3]->GetBinError(b4p);
541  break;
542  }
543  } // for(Int_t b4p=1;b4p<=nBins4p;b4p++)
544  if(TMath::Abs(dCosn1n2n2n1) < 1.e-44)
545  {
546  cout<<Form("labeln1n2n2n1 = %s",labeln1n2n2n1.Data())<<endl;
547  Warning(sMethodName.Data(),"TMath::Abs(dCosn1n2n2n1) < 1.e-44 !!!!");
548  }
549 
550  // Access <<Cos(-n1,n1)>> and <<Cos(-n2,n2)>>:
551  for(Int_t b2p=1;b2p<=nBins2p;b2p++)
552  {
553  if(labeln1n1.EqualTo(fCorrelationsPro[0][1]->GetXaxis()->GetBinLabel(b2p)))
554  {
555  //cout<<labeln1n1.Data()<<endl;
556  dCosn1n1 = fCorrelationsPro[0][1]->GetBinContent(b2p);
557  dCosn1n1Err = fCorrelationsPro[0][1]->GetBinError(b2p);
558  }
559  else if(labeln2n2.EqualTo(fCorrelationsPro[0][1]->GetXaxis()->GetBinLabel(b2p)))
560  {
561  //cout<<labeln2n2.Data()<<endl;
562  dCosn2n2 = fCorrelationsPro[0][1]->GetBinContent(b2p);
563  dCosn2n2Err = fCorrelationsPro[0][1]->GetBinError(b2p);
564  }
565  if(TMath::Abs(dCosn1n1) > 0. && TMath::Abs(dCosn2n2) > 0.){break;} // found 'em both!
566  } // for(Int_t b2p=1;b2p<=nBins2p;b2p++)
567  if(TMath::Abs(dCosn1n1) < 1.e-44)
568  {
569  cout<<Form("labeln1n1 = %s",labeln1n1.Data())<<endl;
570  Warning(sMethodName.Data(),"TMath::Abs(dCosn1n1) < 1.e-44 !!!!");
571  }
572  if(TMath::Abs(dCosn2n2) < 1.e-44)
573  {
574  cout<<Form("labeln2n2 = %s",labeln2n2.Data())<<endl;
575  Warning(sMethodName.Data(),"TMath::Abs(dCosn2n2) < 1.e-44 !!!!");
576  }
577 
578  // Calculate standard candles:
579  dSCn1n2n2n1 = dCosn1n2n2n1-dCosn1n1*dCosn2n2;
580 
581  // Store the final results:
582  fStandardCandlesHist->SetBinContent(b,dSCn1n2n2n1);
583 
584  // Error propagation:
585  if(!fPropagateErrorSC)
586  {
587  fStandardCandlesHist->SetBinError(b,0.);
588  continue;
589  }
590 
591  // Access covariances (multiplied by weight dependent prefactor):
592  Double_t wCovCosn1n2n2n1Cosn1n1 = Covariance(labeln1n2n2n1.Data(),labeln1n1.Data(),fProductsSCPro); // weighted Cov(<Cos(-n1,-n2,n2,n1)>,<Cos(-n1,n1)>)
593  Double_t wCovCosn1n2n2n1Cosn2n2 = Covariance(labeln1n2n2n1.Data(),labeln2n2.Data(),fProductsSCPro); // weighted Cov(<Cos(-n1,-n2,n2,n1)>,<Cos(-n2,n2)>)
594  Double_t wCovCosn1n1Cosn2n2 = Covariance(labeln1n1.Data(),labeln2n2.Data(),fProductsSCPro); // weighted Cov(<Cos(-n1,n1)>,<Cos(-n2,n2)>)
595 
596  // Explicit error propagation:
597  Double_t dSCn1n2n2n1ErrSquared = pow(dCosn1n1,2.)*pow(dCosn2n2Err,2.) + pow(dCosn2n2,2.)*pow(dCosn1n1Err,2.)
598  + pow(dCosn1n2n2n1Err,2.) + 2.*dCosn1n1*dCosn2n2*wCovCosn1n1Cosn2n2
599  - 2.*dCosn1n1*wCovCosn1n2n2n1Cosn2n2 - 2.*dCosn2n2*wCovCosn1n2n2n1Cosn1n1;
600  if(dSCn1n2n2n1ErrSquared > 0.)
601  {
602  dSCn1n2n2n1Err = pow(dSCn1n2n2n1ErrSquared,0.5);
603  } else
604  {
605  Warning(sMethodName.Data(),"dSCn1n2n2n1ErrSquared > 0. is not satisfied for %s !!!!",labeln1n2n2n1.ReplaceAll("Cos","SC").Data());
606  fPropagateError = kFALSE;
607  }
608 
609  // Store the final stat. error:
610  if(fPropagateError)
611  {
612  fStandardCandlesHist->SetBinError(b,dSCn1n2n2n1Err);
613  }
614  } // for(Int_t b=1;b<=nBins;b++)
615 
616  fPropagateError = kTRUE;
617 
618  return;
619 
620 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateStandardCandles()
621 
622 //=======================================================================================================================
623 
625 {
626  // Initialize all arrays for correlations.
627 
628  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
629  {
630  for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
631  {
632  fCorrelationsPro[cs][c] = NULL;
633  }
634  }
635 
636 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForCorrelations()
637 
638 //=======================================================================================================================
639 
641 {
642  // Initialize all arrays for event-by-event cumulants.
643 
644  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
645  {
646  for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
647  {
648  fEbECumulantsPro[cs][c] = NULL;
649  }
650  }
651 
652 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForEbECumulants()
653 
654 //=======================================================================================================================
655 
657 {
658  // Initialize all arrays for differential correlations.
659 
660  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
661  {
662  for(Int_t c=0;c<4;c++) // [1p,2p,3p,4p]
663  {
664  fDiffCorrelationsPro[cs][c] = NULL;
665  fDiffHarmonics[cs][c] = 0;
666  }
667  }
668 
669  // Default values:
670  // Cos, 2p:
671  fDiffHarmonics[1][0] = -2;
672  fDiffHarmonics[1][1] = 2;
673  // Cos, 3p:
674  fDiffHarmonics[2][0] = -3;
675  fDiffHarmonics[2][1] = 1;
676  fDiffHarmonics[2][2] = 2;
677  // Cos, 4p:
678  fDiffHarmonics[3][0] = -2;
679  fDiffHarmonics[3][1] = -2;
680  fDiffHarmonics[3][2] = 2;
681  fDiffHarmonics[3][3] = 2;
682 
683 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForDiffCorrelations()
684 
685 //=======================================================================================================================
686 
688 {
689  // Initialize all arrays for symmetry plane correlations.
690 
691  for(Int_t gc=0;gc<1;gc++) // 'generic correlator': [[0]:(Psi2n,Psi1n),[1]:...] TBI upper boundary will change
692  {
693  for(Int_t n=0;n<2;n++) // 'harmonic n for generic correlator': [[0]:n=1,[1]:n=2,...] TBI upper boundary will change
694  {
695  fSymmetryPlanesPro[gc][n] = NULL;
696  }
697  }
698 
699 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForSymmetryPlanes()
700 
701 //=======================================================================================================================
702 
704 {
705  // Initialize all arrays for nested loops.
706 
707  fCrossCheckDiffCSCOBN[0] = 0; // cos/sin
708  fCrossCheckDiffCSCOBN[1] = 2; // correlator order
709  fCrossCheckDiffCSCOBN[2] = 4; // bin number
710 
711 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForNestedLoops()
712 
713 //=======================================================================================================================
714 
716 {
717  // Initialize all arrays for eta gaps.
718 
719  for(Int_t h=0;h<6;h++) // harmonic
720  {
721  fEtaGapsPro[h] = NULL;
722  }
723 
724 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForEtaGaps()
725 
726 //=======================================================================================================================
727 
729 {
730  // Calculate multi-particle correlations from Q-vector components.
731 
732  // a) Calculate all booked multi-particle correlations;
733  // b) Calculate products needed for QC error propagation;
734  // c) Calculate products needed for SC error propagation.
735 
736  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateCorrelations(AliFlowEventSimple *anEvent)";
737  if(!anEvent){Fatal(sMethodName.Data(),"'anEvent'!?!? You again!!!!");}
738 
739  // a) Calculate all booked multi-particle correlations:
740  Double_t dMultRP = fSelectRandomlyRPs ? fnSelectedRandomlyRPs : anEvent->GetNumberOfRPs(); // TBI shall I promote this variable into data member?
741  if(fSkipSomeIntervals){ dMultRP = dMultRP - fNumberOfSkippedRPParticles; }
742 
743  for(Int_t cs=0;cs<2;cs++) // cos/sin
744  {
745  if(fCalculateOnlyCos && 1==cs){continue;}
746  else if(fCalculateOnlySin && 0==cs){continue;}
747  for(Int_t co=0;co<8;co++) // correlator order (TBI hardwired 8)
748  {
749  if(dMultRP < co+1){break;} // defines min. number of particles in an event for a certain correlator to make sense
750  Int_t nBins = 0;
751  if(fCorrelationsPro[cs][co]){nBins = fCorrelationsPro[cs][co]->GetNbinsX();}
752  else{continue;}
753  for(Int_t b=1;b<=nBins;b++)
754  {
755  TString sBinLabel = fCorrelationsPro[cs][co]->GetXaxis()->GetBinLabel(b);
756  if(sBinLabel.EqualTo("")){break;}
757  Double_t num = CastStringToCorrelation(sBinLabel.Data(),kTRUE);
758  Double_t den = CastStringToCorrelation(sBinLabel.Data(),kFALSE);
759  Double_t weight = den; // TBI: add support for other options for the weight eventually
760  if(den>0.)
761  {
762  fCorrelationsPro[cs][co]->Fill(b-.5,num/den,weight);
763  } else{Warning(sMethodName.Data(),"if(den>0.)");}
764  } // for(Int_t b=1;b<=nBins;b++)
765  } // for(Int_t co=0;co<8;co++) // correlator order (TBI hardwired 8)
766  } // for(Int_t cs=0;cs<=1;cs++) // cos/sin
767 
768  // b) Calculate products needed for QC error propagation:
770 
771  // c) Calculate products needed for SC error propagation:
773 
774 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateCorrelations(AliFlowEventSimple *anEvent)
775 
776 //=======================================================================================================================
777 
779 {
780  // Calculate differential multi-particle correlations from Q-, p- and q-vector components.
781 
782  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateCorrelations(AliFlowEventSimple *anEvent)";
783  if(!anEvent){Fatal(sMethodName.Data(),"'anEvent'!?!? You again!!!!");}
784 
785  Int_t nBins = 0; // TBI promote this to data member?
786  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
787  {
788  if(nBins != 0){break;}
789  for(Int_t co=0;co<4;co++) // [1p,2p,3p,4p]
790  {
791  if(fDiffCorrelationsPro[cs][co] && 0==nBins)
792  {
793  nBins = fDiffCorrelationsPro[cs][co]->GetNbinsX();
794  }
795  } // for(Int_t co=0;co<4;co++) // [1p,2p,3p,4p]
796  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
797 
798  // TBI: The lines below are genuine, most delicious, spaghetti ever... To be reimplemented (one day).
800  {
801  for(Int_t b=1;b<=nBins;b++)
802  {
803  fDiffBinNo = b-1;
804  // <2'>:
805  Double_t num2 = TwoDiff(fDiffHarmonics[1][0],fDiffHarmonics[1][1]).Re();
806  Double_t den2 = TwoDiff(0,0).Re();
807  Double_t w2 = den2; // TBI add support for other options for the weight
808  if(den2>0.){fDiffCorrelationsPro[0][1]->Fill(fDiffCorrelationsPro[0][1]->GetBinCenter(b),num2/den2,w2);}
809  // <3'>:
810  Double_t num3 = ThreeDiff(fDiffHarmonics[2][0],fDiffHarmonics[2][1],fDiffHarmonics[2][2]).Re();
811  Double_t den3 = ThreeDiff(0,0,0).Re();
812  Double_t w3 = den3; // TBI add support for other options for the weight
813  if(den3>0.){fDiffCorrelationsPro[0][2]->Fill(fDiffCorrelationsPro[0][2]->GetBinCenter(b),num3/den3,w3);}
814  // <4'>:
815  Double_t num4 = FourDiff(fDiffHarmonics[3][0],fDiffHarmonics[3][1],fDiffHarmonics[3][2],fDiffHarmonics[3][3]).Re();
816  Double_t den4 = FourDiff(0,0,0,0).Re();
817  Double_t w4 = den4; // TBI add support for other options for the weight
818  if(den4>0.){fDiffCorrelationsPro[0][3]->Fill(fDiffCorrelationsPro[0][3]->GetBinCenter(b),num4/den4,w4);}
819  } // for(Int_t b=1;b<=nBins;b++)
820  }
821  // TBI: The lines below are genuine, most delicious, spaghetti ever... To be reimplemented (one day).
823  {
824  for(Int_t b=1;b<=nBins;b++)
825  {
826  fDiffBinNo = b-1;
827  // <2'>:
828  Double_t num2 = TwoDiff(fDiffHarmonics[1][0],fDiffHarmonics[1][1]).Im();
829  Double_t den2 = TwoDiff(0,0).Re();
830  Double_t w2 = den2; // TBI add support for other options for the weight
831  if(den2>0.){fDiffCorrelationsPro[1][1]->Fill(fDiffCorrelationsPro[1][1]->GetBinCenter(b),num2/den2,w2);}
832  // <3'>:
833  Double_t num3 = ThreeDiff(fDiffHarmonics[2][0],fDiffHarmonics[2][1],fDiffHarmonics[2][2]).Im();
834  Double_t den3 = ThreeDiff(0,0,0).Re();
835  Double_t w3 = den3; // TBI add support for other options for the weight
836  if(den3>0.){fDiffCorrelationsPro[1][2]->Fill(fDiffCorrelationsPro[1][2]->GetBinCenter(b),num3/den3,w3);}
837  // <4'>:
838  Double_t num4 = FourDiff(fDiffHarmonics[3][0],fDiffHarmonics[3][1],fDiffHarmonics[3][2],fDiffHarmonics[3][3]).Im();
839  Double_t den4 = FourDiff(0,0,0,0).Re();
840  Double_t w4 = den4; // TBI add support for other options for the weight
841  if(den4>0.){fDiffCorrelationsPro[1][3]->Fill(fDiffCorrelationsPro[1][3]->GetBinCenter(b),num4/den4,w4);}
842  } // for(Int_t b=1;b<=nBins;b++)
843  }
844 
845 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateDiffCorrelations(AliFlowEventSimple *anEvent)
846 
847 //=======================================================================================================================
848 
850 {
851  // Calculate 2-p correlations with eta gaps.
852 
853  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateEtaGaps(AliFlowEventSimple *anEvent)";
854  if(!anEvent){Fatal(sMethodName.Data(),"'anEvent'? What's wrong with you today...");}
855 
856  TComplex Qa[6][11] = {{TComplex(0.,0.)}}; // -eta [harmonic][eta gap]
857  Double_t Ma[6][11] = {{0.}}; // multiplicity for -eta TBI this shall not depend on harmonic, clearly
858  TComplex Qb[6][11] = {{TComplex(0.,0.)}}; // +eta [harmonic][eta gap]
859  Double_t Mb[6][11] = {{0.}}; // multiplicity for +eta TBI this shall not depend on harmonic, clearly
860  Double_t dEtaGaps[11] = {1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.0};
861 
862  Int_t nTracks = anEvent->NumberOfTracks(); // TBI shall I promote this to data member?
863  Double_t dPhi=0.,dPt=0.,dEta=0.;
864  Double_t wPhi=1.,wPt=1.,wEta=1.;
865  Double_t wToPowerP=1.;
866  for(Int_t t=0;t<nTracks;t++) // loop over all tracks
867  {
868  AliFlowTrackSimple *pTrack = anEvent->GetTrack(t);
869  if(!pTrack){printf("\n pTrack is NULL in MPC::CalculateEtaGaps(AliFlowEventSimple *anEvent) !!!!"); continue;}
870  if(!(pTrack->InRPSelection() || pTrack->InPOISelection())){printf("\n pTrack is neither RP nor POI !!!!"); continue;}
871 
872  if(pTrack->InRPSelection()) // fill Q-vector components only with reference particles
873  {
874  // Access kinematic variables for RP and corresponding weights:
875  dPhi = pTrack->Phi(); // azimuthal angle
876  if(fUseWeights[0][0]){wPhi = Weight(dPhi,"RP","phi");} // corresponding phi weight
877  //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
878  //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
879  dPt = pTrack->Pt();
880  if(fUseWeights[0][1]){wPt = Weight(dPt,"RP","pt");} // corresponding pT weight
881  dEta = pTrack->Eta();
882  if(fUseWeights[0][2]){wEta = Weight(dEta,"RP","eta");} // corresponding eta weight
883  if(fUseWeights[0][0]||fUseWeights[0][1]||fUseWeights[0][2]){wToPowerP = wPhi*wPt*wEta;}
884  // Calculate Qa and Qb vectors:
885  if(dEta<0.) // Qa
886  {
888  {
889  for(Int_t eg=0;eg<11;eg++) // eta gaps
890  {
891  if(dEta<-1.*dEtaGaps[eg]/2.)
892  {
893  Qa[h][eg] += TComplex(wToPowerP*TMath::Cos((h+1)*dPhi),wToPowerP*TMath::Sin((h+1)*dPhi));
894  Ma[h][eg]+=wToPowerP;
895  }
896  } // for(Int_t eg=0;eg<11;eg++) // eta gaps
897  } // for(Int_t h=fLowestHarmonicEtaGaps-1;h<=fHighestHarmonicEtaGaps-1;h++)
898  } // if(dEta<0.)
899  else if(dEta>0.) // Qb
900  {
902  {
903  for(Int_t eg=0;eg<11;eg++) // eta gaps
904  {
905  if(dEta>dEtaGaps[eg]/2.)
906  {
907  Qb[h][eg] += TComplex(wToPowerP*TMath::Cos((h+1)*dPhi),wToPowerP*TMath::Sin((h+1)*dPhi));
908  Mb[h][eg]+=wToPowerP;
909  }
910  } // for(Int_t eg=0;eg<11;eg++) // eta gaps
911  } // for(Int_t h=fLowestHarmonicEtaGaps-1;h<=fHighestHarmonicEtaGaps-1;h++)
912  }
913  } // if(pTrack->InRPSelection())
914  } // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
915 
916  // Calculate 2-p correlations with eta gaps from Qa and Qb vectors:
918  {
919  for(Int_t eg=0;eg<11;eg++) // eta gaps
920  {
921  if(!(Qa[h][eg].Rho()>0. && Qb[h][eg].Rho()>0.)){continue;}
922  if(!(Ma[h][eg]>0. && Mb[h][eg]>0.)){continue;}
923  fEtaGapsPro[h]->Fill(eg+0.5,TComplex(Qa[h][eg]*TComplex::Conjugate(Qb[h][eg])).Re()/(Ma[h][eg]*Mb[h][eg]),Ma[h][eg]*Mb[h][eg]);
924  } // for(Int_t eg=0;eg<11;eg++) // eta gaps
925  } // for(Int_t h=fLowestHarmonicEtaGaps-1;h<=fHighestHarmonicEtaGaps-1;h++)
926 
927 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateEtaGaps(AliFlowEventSimple *anEvent)
928 
929 //=======================================================================================================================
930 
932 {
933  // Cast string of the generic form Cos/Sin(-n_1,-n_2,...,n_{k-1},n_k) in the corresponding correlation value.
934  // If you issue a call to this method with setting numerator = kFALSE, then you are getting back for free
935  // the corresponding denumerator (a.k.a. weight 'number of combinations').
936 
937  // TBI:
938  // a) add protection against cases a la:
939  // string = Cos(-3,-4,5,6,5,6,-3)
940  // method = Six(-3,-4,5,6,5,-3).Re()
941  // b) cross-check with nested loops this method
942 
943  Double_t dValue = 0.; // return value
944 
945  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CastStringToCorrelation(const char *string, Bool_t numerator)";
946 
947  if(!(TString(string).BeginsWith("Cos") || TString(string).BeginsWith("Sin")))
948  {
949  cout<<Form("And the fatal string is... '%s'. Congratulations!!",string)<<endl;
950  Fatal(sMethodName.Data(),"!(TString(string).BeginsWith(...");
951  }
952 
953  Bool_t bRealPart = kTRUE;
954  if(TString(string).BeginsWith("Sin")){bRealPart = kFALSE;}
955 
956  Int_t n[8] = {0,0,0,0,0,0,0,0}; // harmonics, supporting up to 8p correlations
957  UInt_t whichCorr = 0;
958  for(Int_t t=0;t<=TString(string).Length();t++)
959  {
960  if(TString(string[t]).EqualTo(",") || TString(string[t]).EqualTo(")")) // TBI this is just ugly
961  {
962  n[whichCorr] = string[t-1] - '0';
963  if(TString(string[t-2]).EqualTo("-")){n[whichCorr] = -1*n[whichCorr];}
964  if(!(TString(string[t-2]).EqualTo("-")
965  || TString(string[t-2]).EqualTo(",")
966  || TString(string[t-2]).EqualTo("("))) // TBI relax this eventually to allow two-digits harmonics
967  {
968  cout<<Form("And the fatal string is... '%s'. Congratulations!!",string)<<endl;
969  Fatal(sMethodName.Data(),"!(TString(string[t-2]).EqualTo(...");
970  }
971  whichCorr++;
972  if(whichCorr>=9){Fatal(sMethodName.Data(),"whichCorr>=9");} // not supporting corr. beyond 8p
973  } // if(TString(string[t]).EqualTo(",") || TString(string[t]).EqualTo(")")) // TBI this is just ugly
974  } // for(UInt_t t=0;t<=TString(string).Length();t++)
975 
976  switch(whichCorr)
977  {
978  case 1:
979  if(!numerator){dValue = One(0).Re();}
980  else if(bRealPart){dValue = One(n[0]).Re();}
981  else{dValue = One(n[0]).Im();}
982  break;
983 
984  case 2:
985  if(!numerator){dValue = Two(0,0).Re();}
986  else if(bRealPart){dValue = Two(n[0],n[1]).Re();}
987  else{dValue = Two(n[0],n[1]).Im();}
988  break;
989 
990  case 3:
991  if(!numerator){dValue = Three(0,0,0).Re();}
992  else if(bRealPart){dValue = Three(n[0],n[1],n[2]).Re();}
993  else{dValue = Three(n[0],n[1],n[2]).Im();}
994  break;
995 
996  case 4:
997  if(!numerator){dValue = Four(0,0,0,0).Re();}
998  else if(bRealPart){dValue = Four(n[0],n[1],n[2],n[3]).Re();}
999  else{dValue = Four(n[0],n[1],n[2],n[3]).Im();}
1000  break;
1001 
1002  case 5:
1003  if(!numerator){dValue = Five(0,0,0,0,0).Re();}
1004  else if(bRealPart){dValue = Five(n[0],n[1],n[2],n[3],n[4]).Re();}
1005  else{dValue = Five(n[0],n[1],n[2],n[3],n[4]).Im();}
1006  break;
1007 
1008  case 6:
1009  if(!numerator){dValue = Six(0,0,0,0,0,0).Re();}
1010  else if(bRealPart){dValue = Six(n[0],n[1],n[2],n[3],n[4],n[5]).Re();}
1011  else{dValue = Six(n[0],n[1],n[2],n[3],n[4],n[5]).Im();}
1012  break;
1013 
1014  case 7:
1015  if(!numerator){dValue = Seven(0,0,0,0,0,0,0).Re();}
1016  else if(bRealPart){dValue = Seven(n[0],n[1],n[2],n[3],n[4],n[5],n[6]).Re();}
1017  else{dValue = Seven(n[0],n[1],n[2],n[3],n[4],n[5],n[6]).Im();}
1018  break;
1019 
1020  case 8:
1021  if(!numerator){dValue = Eight(0,0,0,0,0,0,0,0).Re();}
1022  else if(bRealPart){dValue = Eight(n[0],n[1],n[2],n[3],n[4],n[5],n[6],n[7]).Re();}
1023  else{dValue = Eight(n[0],n[1],n[2],n[3],n[4],n[5],n[6],n[7]).Im();}
1024  break;
1025 
1026  default:
1027  cout<<Form("And the fatal 'whichCorr' value is... %d. Congratulations!!",whichCorr)<<endl;
1028  Fatal(sMethodName.Data(),"switch(whichCorr)");
1029  } // switch(whichCorr)
1030 
1031  return dValue;
1032 
1033 } // Double_t AliFlowAnalysisWithMultiparticleCorrelations::CastStringToCorrelation(const char *string, Bool_t numerator)
1034 
1035 //=======================================================================================================================
1036 
1038 {
1039  // Calculate products of multi-particle correlations (needed for error propagation).
1040 
1041  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateProductsOfCorrelations(AliFlowEventSimple *anEvent, TProfile2D *profile2D)";
1042  if(!anEvent){Fatal(sMethodName.Data(),"Sorry, 'anEvent' is on holidays.");}
1043  if(!profile2D){Fatal(sMethodName.Data(),"Sorry, 'profile2D' is on holidays.");}
1044 
1045  Int_t nBins = profile2D->GetXaxis()->GetNbins();
1046  for(Int_t bx=2;bx<=nBins;bx++)
1047  {
1048  for(Int_t by=1;by<bx;by++)
1049  {
1050  const char *binLabelX = profile2D->GetXaxis()->GetBinLabel(bx);
1051  const char *binLabelY = profile2D->GetYaxis()->GetBinLabel(by);
1052  Double_t numX = this->CastStringToCorrelation(binLabelX,kTRUE); // numerator
1053  Double_t denX = this->CastStringToCorrelation(binLabelX,kFALSE); // denominator
1054  Double_t wX = denX; // weight TBI add support for other options
1055  Double_t numY = this->CastStringToCorrelation(binLabelY,kTRUE); // numerator
1056  Double_t denY = this->CastStringToCorrelation(binLabelY,kFALSE); // denominator
1057  Double_t wY = denY; // weight TBI add support for other options
1058  if(TMath::Abs(denX) > 0. && TMath::Abs(denY) > 0.)
1059  {
1060  profile2D->Fill(bx-0.5,by-0.5,(numX/denX)*(numY/denY),wX*wY);
1061  } else
1062  {
1063  cout<<endl;
1064  cout<<"Cannot calculate product for:"<<endl;
1065  cout<<Form("binLabelX = %s",binLabelX)<<endl;
1066  cout<<Form("binLabelY = %s",binLabelY)<<endl;
1067  cout<<Form("anEvent->GetNumberOfRPs() = %d",anEvent->GetNumberOfRPs())<<endl;
1068  cout<<Form("fNumberOfSkippedRPParticles = %d",fNumberOfSkippedRPParticles)<<endl;
1069  Fatal(sMethodName.Data(),"if(TMath::Abs(denX) > 0. && TMath::Abs(denY) > 0.)");
1070  } // else
1071  } // for(Int_t by=1;by<bx;by++)
1072  } // for(Int_t bx=2;bx<=nBins;bx++)
1073 
1074 } // void CalculateProductsOfCorrelations(AliFlowEventSimple *anEvent, TProfile2D *profile2D)
1075 
1076 //=======================================================================================================================
1077 
1079 {
1080  // Calculate e-b-e cumulants from Q-vector components.
1081 
1082  // TBI this mathod is far (very far, in fact) from being finalized :'(
1083 
1084  // a) Calculate and store e-b-e cumulants.
1085 
1086  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CalculateEbECumulants(AliFlowEventSimple *anEvent)";
1087  if(!anEvent){Fatal(sMethodName.Data(),"'anEvent'!?!? You again!!!!");}
1088 
1089  // a) Calculate and store e-b-e cumulants:
1090  Double_t dMultRP = anEvent->GetNumberOfRPs(); // TBI shall I promote this variable into data member?
1091 
1092  if(fSkipSomeIntervals){ dMultRP = dMultRP - fNumberOfSkippedRPParticles; } // TBI tmp gym
1093 
1094  Int_t binNo[8]; for(Int_t c=0;c<8;c++){binNo[c]=1;}
1095  // 1-p:
1096  for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++)
1097  {
1098  if(fSkipZeroHarmonics && 0==n1){continue;}
1099  if(fCalculateAll)
1100  {
1101  TComplex oneN = One(n1); // numerator
1102  Double_t oneD = One(0).Re(); // denominator
1103  Double_t oneW = oneD; // weight TBI add other possibilities here for the weight
1104  if(oneD>0. && dMultRP>=1)
1105  {
1106  fEbECumulantsPro[0][0]->Fill(binNo[0]-.5,oneN.Re()/oneD,oneW);
1107  fEbECumulantsPro[1][0]->Fill(binNo[0]++-.5,oneN.Im()/oneD,oneW);
1108  } else {Warning(sMethodName.Data(),"if(oneD>0. && dMultRP>=1) ");}
1109  }
1110  if(1==fDontGoBeyond){continue;}
1111  // 2-p:
1112  for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
1113  {
1114  if(fSkipZeroHarmonics && 0==n2){continue;}
1115  if(fCalculateAll
1116  || (fCalculateIsotropic && 0==n1+n2)
1117  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2))
1118  || (fCalculateSameIsotropic && 0==n1+n2 && TMath::Abs(n1)==TMath::Abs(n2)))
1119  {
1120  Double_t cumulants2pCos = Two(n1,n2).Re()/Two(0,0).Re()
1121  - (One(n1).Re()/One(0).Re())*(One(n2).Re()/One(0).Re())
1122  + (One(n1).Im()/One(0).Re())*(One(n2).Im()/One(0).Re());
1123 
1124  Double_t cumulants2pSin = Two(n1,n2).Im()/Two(0,0).Re()
1125  - (One(n1).Re()/One(0).Re())*(One(n2).Im()/One(0).Re())
1126  - (One(n2).Re()/One(0).Re())*(One(n1).Im()/One(0).Re());
1127 
1128  if(/*twoD>0. &&*/ dMultRP>=2)
1129  {
1130  fEbECumulantsPro[0][1]->Fill(binNo[1]-.5,cumulants2pCos,1.);;
1131  fEbECumulantsPro[1][1]->Fill(binNo[1]++-.5,cumulants2pSin,1.);;
1132  } else {Warning(sMethodName.Data(),"/*twoD>0. &&*/ dMultRP>=2");}
1133  }
1134  if(2==fDontGoBeyond){continue;}
1135 
1136  /*
1137 
1138  // 3-p:
1139  for(Int_t n3=n2;n3<=fMaxHarmonic;n3++)
1140  {
1141  if(fSkipZeroHarmonics && 0==n3){continue;}
1142  if(fCalculateAll
1143  || (fCalculateIsotropic && 0==n1+n2+n3)
1144  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3))
1145  || (fCalculateSameIsotropic && 0==n1+n2+n3 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)))
1146  {
1147  TComplex threeN = Three(n1,n2,n3); // numerator
1148  Double_t threeD = Three(0,0,0).Re(); // denominator
1149  Double_t threeW = threeD; // weight TBI add other possibilities here for the weight
1150  if(threeD>0. && dMultRP>=3)
1151  {
1152  fEbECumulantsPro[0][2]->Fill(binNo[2]-.5,threeN.Re()/threeD,threeW);
1153  fEbECumulantsPro[1][2]->Fill(binNo[2]++-.5,threeN.Im()/threeD,threeW);
1154  } else {Warning(sMethodName.Data(),"threeD>0. && dMultRP>=3");}
1155  }
1156  if(3==fDontGoBeyond){continue;}
1157  // 4-p:
1158  for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)
1159  {
1160  if(fSkipZeroHarmonics && 0==n4){continue;}
1161  if(fCalculateAll
1162  || (fCalculateIsotropic && 0==n1+n2+n3+n4)
1163  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
1164  && TMath::Abs(n1)==TMath::Abs(n4))
1165  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4 && TMath::Abs(n1)==TMath::Abs(n2)
1166  && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)))
1167  {
1168  TComplex fourN = Four(n1,n2,n3,n4); // numerator
1169  Double_t fourD = Four(0,0,0,0).Re(); // denominator
1170  Double_t fourW = fourD; // weight TBI add other possibilities here for the weight
1171  if(fourD>0. && dMultRP>=4)
1172  {
1173  fEbECumulantsPro[0][3]->Fill(binNo[3]-.5,fourN.Re()/fourD,fourW);
1174  fEbECumulantsPro[1][3]->Fill(binNo[3]++-.5,fourN.Im()/fourD,fourW);
1175  } else {Warning(sMethodName.Data(),"fourD>0. && dMultRP>=4");}
1176  }
1177  if(4==fDontGoBeyond){continue;}
1178  // 5-p:
1179  for(Int_t n5=n4;n5<=fMaxHarmonic;n5++)
1180  {
1181  if(fSkipZeroHarmonics && 0==n5){continue;}
1182  if(fCalculateAll
1183  || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5)
1184  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
1185  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5))
1186  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5 && TMath::Abs(n1)==TMath::Abs(n2)
1187  && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5)))
1188  {
1189  TComplex fiveN = Five(n1,n2,n3,n4,n5); // numerator
1190  Double_t fiveD = Five(0,0,0,0,0).Re(); // denominator
1191  Double_t fiveW = fiveD; // weight TBI add other possibilities here for the weight
1192  if(fiveD>0. && dMultRP>=5)
1193  {
1194  fEbECumulantsPro[0][4]->Fill(binNo[4]-.5,fiveN.Re()/fiveD,fiveW);
1195  fEbECumulantsPro[1][4]->Fill(binNo[4]++-.5,fiveN.Im()/fiveD,fiveW);
1196  } else {Warning(sMethodName.Data(),"fiveD>0. && dMultRP>=5");}
1197  }
1198  if(5==fDontGoBeyond){continue;}
1199  // 6-p:
1200  for(Int_t n6=n5;n6<=fMaxHarmonic;n6++)
1201  {
1202  if(fSkipZeroHarmonics && 0==n6){continue;}
1203  if(fCalculateAll
1204  || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6)
1205  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
1206  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6))
1207  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
1208  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6)))
1209  {
1210  TComplex sixN = Six(n1,n2,n3,n4,n5,n6); // numerator
1211  Double_t sixD = Six(0,0,0,0,0,0).Re(); // denominator
1212  Double_t sixW = sixD; // weight TBI add other possibilities here for the weight
1213  if(sixD>0. && dMultRP>=6)
1214  {
1215  fEbECumulantsPro[0][5]->Fill(binNo[5]-.5,sixN.Re()/sixD,sixW);
1216  fEbECumulantsPro[1][5]->Fill(binNo[5]++-.5,sixN.Im()/sixD,sixW);
1217  } else {Warning(sMethodName.Data(),"sixD>0. && dMultRP>=6");}
1218  }
1219  if(6==fDontGoBeyond){continue;}
1220  // 7-p:
1221  for(Int_t n7=n6;n7<=fMaxHarmonic;n7++)
1222  {
1223  if(fSkipZeroHarmonics && 0==n7){continue;}
1224  if(fCalculateAll
1225  || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7)
1226  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
1227  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7))
1228  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
1229  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6)
1230  && TMath::Abs(n1)==TMath::Abs(n7)))
1231  {
1232  TComplex sevenN = Seven(n1,n2,n3,n4,n5,n6,n7); // numerator
1233  Double_t sevenD = Seven(0,0,0,0,0,0,0).Re(); // denominator
1234  Double_t sevenW = sevenD; // weight TBI add other possibilities here for the weight
1235  if(sevenD>0. && dMultRP>=7)
1236  {
1237  fEbECumulantsPro[0][6]->Fill(binNo[6]-.5,sevenN.Re()/sevenD,sevenW);
1238  fEbECumulantsPro[1][6]->Fill(binNo[6]++-.5,sevenN.Im()/sevenD,sevenW);
1239  } else {Warning(sMethodName.Data(),"sevenD>0. && dMultRP>=7");}
1240  }
1241  if(7==fDontGoBeyond){continue;}
1242  // 8-p:
1243  for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
1244  {
1245  if(fSkipZeroHarmonics && 0==n8){continue;}
1246  if(fCalculateAll
1247  || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8)
1248  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
1249  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7)
1250  && TMath::Abs(n1)==TMath::Abs(n8))
1251  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
1252  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6)
1253  && TMath::Abs(n1)==TMath::Abs(n7) && TMath::Abs(n1)==TMath::Abs(n8)))
1254  {
1255  TComplex eightN = Eight(n1,n2,n3,n4,n5,n6,n7,n8); // numerator
1256  Double_t eightD = Eight(0,0,0,0,0,0,0,0).Re(); // denominator
1257  Double_t eightW = eightD; // weight TBI add other possibilities here for the weight
1258  if(eightD>0. && dMultRP>=8)
1259  {
1260  fEbECumulantsPro[0][7]->Fill(binNo[7]-.5,eightN.Re()/eightD,eightW);
1261  fEbECumulantsPro[1][7]->Fill(binNo[7]++-.5,eightN.Im()/eightD,eightW);
1262  }
1263  }
1264  } // for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
1265  } // for(Int_t n7=n6;n7<=fMaxHarmonic;n7++)
1266  } // for(Int_t n6=n5;n6<=fMaxHarmonic;n6++)
1267  } // for(Int_t n5=n4;n5<=fMaxHarmonic;n5++)
1268  } // for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)
1269  } // for(Int_t n3=n2;n3<=fMaxHarmonic;n3++)
1270 
1271  */
1272 
1273  } // for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
1274  } // for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++)
1275 
1276 } // void AliFlowAnalysisWithMultiparticleCorrelations::CalculateEbECumulants(AliFlowEventSimple *anEvent)
1277 
1278 //=======================================================================================================================
1279 
1281 {
1282  // TBI
1283 
1284  if(!pTrack){exit(0);} // TBI
1285 
1286  Double_t dPhi = pTrack->Phi();
1287  Double_t dPt = pTrack->Pt();
1288  Double_t dEta = pTrack->Eta();
1289  Double_t dPhiPtEta[3] = {dPhi,dPt,dEta};
1290 
1291  // Skip some intervals: TBI promote eventually to AFTC class
1292  Bool_t bPasses = kTRUE;
1293  for(Int_t ppe=0;ppe<3;ppe++)
1294  {
1295  if(!bPasses){break;} // found one kinematic window which shall be skipped
1296  for(Int_t b=0;b<10;b+=2)
1297  {
1298  if(-44==(Int_t)fSkip[ppe][b]){continue;}
1299  if(dPhiPtEta[ppe]>=fSkip[ppe][b] && dPhiPtEta[ppe]<fSkip[ppe][b+1])
1300  {
1301  bPasses = kFALSE;
1302  break;
1303  } // TBI this is a clear bug when this setter is used multiple times...
1304  } // for(Int_t b=0;b<10;b++)
1305  } // for(Int_t ppe=0;ppe<3;ppe++)
1306 
1307  return bPasses;
1308 
1309 } // Bool_t AliFlowAnalysisWithMultiparticleCorrelations::TrackIsInSpecifiedIntervals(AliFlowTrackSimple *pTrack)
1310 
1311 //=======================================================================================================================
1312 
1314 {
1315  // Cross-check results for multi-particle correlations with nested loops.
1316 
1317  // TBI add few comments here, there and over there
1318  // TBI this method is rather messy :'(
1319 
1320  Int_t h1 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(2);
1321  Int_t h2 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(3);
1322  Int_t h3 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(4);
1323  Int_t h4 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(5);
1324  Int_t h5 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(6);
1325  Int_t h6 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(7);
1326  Int_t h7 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(8);
1327  Int_t h8 = (Int_t)fNestedLoopsFlagsPro->GetBinContent(9);
1328 
1329  this->ResetQvector();
1330  this->FillQvector(anEvent);
1331 
1332  if(TMath::Abs(One(0).Re())>0.)
1333  {
1334  fNestedLoopsResultsCosPro->Fill(1.5,One(h1).Re()/One(0).Re(),One(0).Re());
1335  fNestedLoopsResultsSinPro->Fill(1.5,One(h1).Im()/One(0).Re(),One(0).Re());
1336  }
1337  if(TMath::Abs(Two(0,0).Re())>0.)
1338  {
1339  fNestedLoopsResultsCosPro->Fill(3.5,Two(h1,h2).Re()/Two(0,0).Re(),Two(0,0).Re());
1340  fNestedLoopsResultsSinPro->Fill(3.5,Two(h1,h2).Im()/Two(0,0).Re(),Two(0,0).Re());
1341  }
1342  if(TMath::Abs(Three(0,0,0).Re())>0.)
1343  {
1344  fNestedLoopsResultsCosPro->Fill(5.5,Three(h1,h2,h3).Re()/Three(0,0,0).Re(),Three(0,0,0).Re());
1345  fNestedLoopsResultsSinPro->Fill(5.5,Three(h1,h2,h3).Im()/Three(0,0,0).Re(),Three(0,0,0).Re());
1346  }
1347  if(TMath::Abs(Four(0,0,0,0).Re())>0.)
1348  {
1349  fNestedLoopsResultsCosPro->Fill(7.5,Four(h1,h2,h3,h4).Re()/Four(0,0,0,0).Re(),Four(0,0,0,0).Re());
1350  fNestedLoopsResultsSinPro->Fill(7.5,Four(h1,h2,h3,h4).Im()/Four(0,0,0,0).Re(),Four(0,0,0,0).Re());
1351  }
1352  if(TMath::Abs(Five(0,0,0,0,0).Re())>0.)
1353  {
1354  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());
1355  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());
1356  }
1357  if(TMath::Abs(Six(0,0,0,0,0,0).Re())>0.)
1358  {
1359  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());
1360  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());
1361  }
1362  if(TMath::Abs(Seven(0,0,0,0,0,0,0).Re())>0.)
1363  {
1364  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());
1365  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());
1366  }
1367  if(TMath::Abs(Eight(0,0,0,0,0,0,0,0).Re())>0.)
1368  {
1369  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());
1370  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());
1371  }
1372 
1373  Int_t nPrim = anEvent->NumberOfTracks();
1374  Double_t dMultRP = anEvent->GetNumberOfRPs(); // TBI shall I promote this variable into data member?
1375  AliFlowTrackSimple *aftsTrack = NULL;
1376  Double_t dPhi1=0.,dPhi2=0.,dPhi3=0.,dPhi4=0.,dPhi5=0.,dPhi6=0.,dPhi7=0.,dPhi8=0.;
1377  Double_t wPhi1=1.,wPhi2=1.,wPhi3=1.,wPhi4=1.,wPhi5=1.,wPhi6=1.,wPhi7=1.,wPhi8=1.;
1378 
1379  // 1-particle stuff: TBI
1380  if(dMultRP>=1)
1381  {
1382  for(Int_t i1=0;i1<nPrim;i1++)
1383  {
1384  aftsTrack = anEvent->GetTrack(i1);
1385  if(!(aftsTrack->InRPSelection())){continue;}
1386  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1387  dPhi1 = aftsTrack->Phi();
1388  if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1389  // Fill:
1390  fNestedLoopsResultsCosPro->Fill(0.5,TMath::Cos(h1*dPhi1),wPhi1);
1391  fNestedLoopsResultsSinPro->Fill(0.5,TMath::Sin(h1*dPhi1),wPhi1);
1392  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1393  } // end of if(nPrim>=1)
1394 
1395  // 2-particle correlations:
1396  if(dMultRP>=2)
1397  {
1398  for(Int_t i1=0;i1<nPrim;i1++)
1399  {
1400  aftsTrack = anEvent->GetTrack(i1);
1401  if(!(aftsTrack->InRPSelection())){continue;}
1402  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1403  dPhi1 = aftsTrack->Phi();
1404  if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1405  for(Int_t i2=0;i2<nPrim;i2++)
1406  {
1407  if(i2==i1){continue;}
1408  aftsTrack = anEvent->GetTrack(i2);
1409  if(!(aftsTrack->InRPSelection())){continue;}
1410  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1411  dPhi2 = aftsTrack->Phi();
1412  if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1413  // Fill:
1414  fNestedLoopsResultsCosPro->Fill(2.5,TMath::Cos(h1*dPhi1+h2*dPhi2),wPhi1*wPhi2);
1415  fNestedLoopsResultsSinPro->Fill(2.5,TMath::Sin(h1*dPhi1+h2*dPhi2),wPhi1*wPhi2);
1416  } // end of for(Int_t i2=0;i2<nPrim;i2++)
1417  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1418  } // end of if(nPrim>=2)
1419 
1420  // 3-particle correlations:
1421  if(dMultRP>=3)
1422  {
1423  for(Int_t i1=0;i1<nPrim;i1++)
1424  {
1425  aftsTrack=anEvent->GetTrack(i1);
1426  if(!(aftsTrack->InRPSelection())){continue;}
1427  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1428  dPhi1=aftsTrack->Phi();
1429  if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1430  for(Int_t i2=0;i2<nPrim;i2++)
1431  {
1432  if(i2==i1){continue;}
1433  aftsTrack=anEvent->GetTrack(i2);
1434  if(!(aftsTrack->InRPSelection())){continue;}
1435  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1436  dPhi2=aftsTrack->Phi();
1437  if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1438  for(Int_t i3=0;i3<nPrim;i3++)
1439  {
1440  if(i3==i1||i3==i2){continue;}
1441  aftsTrack=anEvent->GetTrack(i3);
1442  if(!(aftsTrack->InRPSelection())){continue;}
1443  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1444  dPhi3=aftsTrack->Phi();
1445  if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1446  // Fill:
1447  fNestedLoopsResultsCosPro->Fill(4.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3),wPhi1*wPhi2*wPhi3);
1448  fNestedLoopsResultsSinPro->Fill(4.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3),wPhi1*wPhi2*wPhi3);
1449  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1450  } // end of for(Int_t i2=0;i2<nPrim;i2++)
1451  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1452  } // end of if(nPrim>=3)
1453 
1454  // 4-particle correlations:
1455  if(dMultRP>=4)
1456  {
1457  for(Int_t i1=0;i1<nPrim;i1++)
1458  {
1459  aftsTrack=anEvent->GetTrack(i1);
1460  if(!(aftsTrack->InRPSelection())){continue;}
1461  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1462  dPhi1=aftsTrack->Phi();
1463  if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1464  for(Int_t i2=0;i2<nPrim;i2++)
1465  {
1466  if(i2==i1){continue;}
1467  aftsTrack=anEvent->GetTrack(i2);
1468  if(!(aftsTrack->InRPSelection())){continue;}
1469  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1470  dPhi2=aftsTrack->Phi();
1471  if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1472  for(Int_t i3=0;i3<nPrim;i3++)
1473  {
1474  if(i3==i1||i3==i2){continue;}
1475  aftsTrack=anEvent->GetTrack(i3);
1476  if(!(aftsTrack->InRPSelection())){continue;}
1477  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1478  dPhi3=aftsTrack->Phi();
1479  if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1480  for(Int_t i4=0;i4<nPrim;i4++)
1481  {
1482  if(i4==i1||i4==i2||i4==i3){continue;}
1483  aftsTrack=anEvent->GetTrack(i4);
1484  if(!(aftsTrack->InRPSelection())){continue;}
1485  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1486  dPhi4=aftsTrack->Phi();
1487  if(fUseWeights[0][0]){wPhi4 = Weight(dPhi4,"RP","phi");}
1488  // Fill:
1489  fNestedLoopsResultsCosPro->Fill(6.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4),wPhi1*wPhi2*wPhi3*wPhi4);
1490  fNestedLoopsResultsSinPro->Fill(6.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4),wPhi1*wPhi2*wPhi3*wPhi4);
1491  } // end of for(Int_t i4=0;i4<nPrim;i4++)
1492  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1493  } // end of for(Int_t i2=0;i2<nPrim;i2++)
1494  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1495  } // end of if(nPrim>=)
1496 
1497  // 5-particle correlations:
1498  if(dMultRP>=5)
1499  {
1500  for(Int_t i1=0;i1<nPrim;i1++)
1501  {
1502  aftsTrack=anEvent->GetTrack(i1);
1503  if(!(aftsTrack->InRPSelection())){continue;}
1504  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1505  dPhi1=aftsTrack->Phi();
1506  if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1507  for(Int_t i2=0;i2<nPrim;i2++)
1508  {
1509  if(i2==i1){continue;}
1510  aftsTrack=anEvent->GetTrack(i2);
1511  if(!(aftsTrack->InRPSelection())){continue;}
1512  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1513  dPhi2=aftsTrack->Phi();
1514  if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1515  for(Int_t i3=0;i3<nPrim;i3++)
1516  {
1517  if(i3==i1||i3==i2){continue;}
1518  aftsTrack=anEvent->GetTrack(i3);
1519  if(!(aftsTrack->InRPSelection())){continue;}
1520  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1521  dPhi3=aftsTrack->Phi();
1522  if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1523  for(Int_t i4=0;i4<nPrim;i4++)
1524  {
1525  if(i4==i1||i4==i2||i4==i3){continue;}
1526  aftsTrack=anEvent->GetTrack(i4);
1527  if(!(aftsTrack->InRPSelection())){continue;}
1528  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1529  dPhi4=aftsTrack->Phi();
1530  if(fUseWeights[0][0]){wPhi4 = Weight(dPhi4,"RP","phi");}
1531  for(Int_t i5=0;i5<nPrim;i5++)
1532  {
1533  if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1534  aftsTrack=anEvent->GetTrack(i5);
1535  if(!(aftsTrack->InRPSelection())){continue;}
1536  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1537  dPhi5=aftsTrack->Phi();
1538  if(fUseWeights[0][0]){wPhi5 = Weight(dPhi5,"RP","phi");}
1539  // Fill:
1540  fNestedLoopsResultsCosPro->Fill(8.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5);
1541  fNestedLoopsResultsSinPro->Fill(8.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5);
1542  } // end of for(Int_t i5=0;i5<nPrim;i5++)
1543  } // end of for(Int_t i4=0;i4<nPrim;i4++)
1544  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1545  } // end of for(Int_t i2=0;i2<nPrim;i2++)
1546  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1547  } // end of if(nPrim>=5)
1548 
1549  // 6-particle correlations:
1550  if(dMultRP>=6)
1551  {
1552  for(Int_t i1=0;i1<nPrim;i1++)
1553  {
1554  aftsTrack=anEvent->GetTrack(i1);
1555  if(!(aftsTrack->InRPSelection())){continue;}
1556  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1557  dPhi1=aftsTrack->Phi();
1558  if(fUseWeights[0][0]){wPhi1 = Weight(dPhi1,"RP","phi");}
1559  for(Int_t i2=0;i2<nPrim;i2++)
1560  {
1561  if(i2==i1){continue;}
1562  aftsTrack=anEvent->GetTrack(i2);
1563  if(!(aftsTrack->InRPSelection())){continue;}
1564  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1565  dPhi2=aftsTrack->Phi();
1566  if(fUseWeights[0][0]){wPhi2 = Weight(dPhi2,"RP","phi");}
1567  for(Int_t i3=0;i3<nPrim;i3++)
1568  {
1569  if(i3==i1||i3==i2){continue;}
1570  aftsTrack=anEvent->GetTrack(i3);
1571  if(!(aftsTrack->InRPSelection())){continue;}
1572  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1573  dPhi3=aftsTrack->Phi();
1574  if(fUseWeights[0][0]){wPhi3 = Weight(dPhi3,"RP","phi");}
1575  for(Int_t i4=0;i4<nPrim;i4++)
1576  {
1577  if(i4==i1||i4==i2||i4==i3){continue;}
1578  aftsTrack=anEvent->GetTrack(i4);
1579  if(!(aftsTrack->InRPSelection())){continue;}
1580  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1581  dPhi4=aftsTrack->Phi();
1582  if(fUseWeights[0][0]){wPhi4 = Weight(dPhi4,"RP","phi");}
1583  for(Int_t i5=0;i5<nPrim;i5++)
1584  {
1585  if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1586  aftsTrack=anEvent->GetTrack(i5);
1587  if(!(aftsTrack->InRPSelection())){continue;}
1588  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1589  dPhi5=aftsTrack->Phi();
1590  if(fUseWeights[0][0]){wPhi5=Weight(dPhi5,"RP","phi");}
1591  for(Int_t i6=0;i6<nPrim;i6++)
1592  {
1593  if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5){continue;}
1594  aftsTrack=anEvent->GetTrack(i6);
1595  if(!(aftsTrack->InRPSelection())){continue;}
1596  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1597  dPhi6=aftsTrack->Phi();
1598  if(fUseWeights[0][0]){wPhi6=Weight(dPhi6,"RP","phi");}
1599  // Fill:
1600  fNestedLoopsResultsCosPro->Fill(10.5,TMath::Cos(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6);
1601  fNestedLoopsResultsSinPro->Fill(10.5,TMath::Sin(h1*dPhi1+h2*dPhi2+h3*dPhi3+h4*dPhi4+h5*dPhi5+h6*dPhi6),wPhi1*wPhi2*wPhi3*wPhi4*wPhi5*wPhi6);
1602  } // end of for(Int_t i6=0;i6<nPrim;i6++)
1603  } // end of for(Int_t i5=0;i5<nPrim;i5++)
1604  } // end of for(Int_t i4=0;i4<nPrim;i4++)
1605  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1606  } // end of for(Int_t i2=0;i2<nPrim;i2++)
1607  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1608  } // end of if(nPrim>=6)
1609 
1610  // 7-particle correlations:
1611  if(dMultRP>=7)
1612  {
1613  for(Int_t i1=0;i1<nPrim;i1++)
1614  {
1615  aftsTrack=anEvent->GetTrack(i1);
1616  if(!(aftsTrack->InRPSelection())){continue;}
1617  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1618  dPhi1=aftsTrack->Phi();
1619  if(fUseWeights[0][0]){wPhi1=Weight(dPhi1,"RP","phi");}
1620  for(Int_t i2=0;i2<nPrim;i2++)
1621  {
1622  if(i2==i1){continue;}
1623  aftsTrack=anEvent->GetTrack(i2);
1624  if(!(aftsTrack->InRPSelection())){continue;}
1625  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1626  dPhi2=aftsTrack->Phi();
1627  if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1628  for(Int_t i3=0;i3<nPrim;i3++)
1629  {
1630  if(i3==i1||i3==i2){continue;}
1631  aftsTrack=anEvent->GetTrack(i3);
1632  if(!(aftsTrack->InRPSelection())){continue;}
1633  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1634  dPhi3=aftsTrack->Phi();
1635  if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1636  for(Int_t i4=0;i4<nPrim;i4++)
1637  {
1638  if(i4==i1||i4==i2||i4==i3){continue;}
1639  aftsTrack=anEvent->GetTrack(i4);
1640  if(!(aftsTrack->InRPSelection())){continue;}
1641  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1642  dPhi4=aftsTrack->Phi();
1643  if(fUseWeights[0][0]){wPhi4=Weight(dPhi4,"RP","phi");}
1644  for(Int_t i5=0;i5<nPrim;i5++)
1645  {
1646  if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1647  aftsTrack=anEvent->GetTrack(i5);
1648  if(!(aftsTrack->InRPSelection())){continue;}
1649  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1650  dPhi5=aftsTrack->Phi();
1651  if(fUseWeights[0][0]){wPhi5=Weight(dPhi5,"RP","phi");}
1652  for(Int_t i6=0;i6<nPrim;i6++)
1653  {
1654  if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5){continue;}
1655  aftsTrack=anEvent->GetTrack(i6);
1656  if(!(aftsTrack->InRPSelection())){continue;}
1657  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1658  dPhi6=aftsTrack->Phi();
1659  if(fUseWeights[0][0]){wPhi6=Weight(dPhi6,"RP","phi");}
1660  for(Int_t i7=0;i7<nPrim;i7++)
1661  {
1662  if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6){continue;}
1663  aftsTrack=anEvent->GetTrack(i7);
1664  if(!(aftsTrack->InRPSelection())){continue;}
1665  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1666  dPhi7=aftsTrack->Phi();
1667  if(fUseWeights[0][0]){wPhi7=Weight(dPhi7,"RP","phi");}
1668  // Fill:
1669  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);
1670  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);
1671  } // end of for(Int_t i7=0;i7<nPrim;i7++)
1672  } // end of for(Int_t i6=0;i6<nPrim;i6++)
1673  } // end of for(Int_t i5=0;i5<nPrim;i5++)
1674  } // end of for(Int_t i4=0;i4<nPrim;i4++)
1675  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1676  } // end of for(Int_t i2=0;i2<nPrim;i2++)
1677  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1678  } // end of if(nPrim>=7)
1679 
1680  // 8-particle correlations:
1681  if(dMultRP>=8)
1682  {
1683  for(Int_t i1=0;i1<nPrim;i1++)
1684  {
1685  aftsTrack=anEvent->GetTrack(i1);
1686  if(!(aftsTrack->InRPSelection())){continue;}
1687  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1688  dPhi1=aftsTrack->Phi();
1689  if(fUseWeights[0][0]){wPhi1=Weight(dPhi1,"RP","phi");}
1690  for(Int_t i2=0;i2<nPrim;i2++)
1691  {
1692  if(i2==i1){continue;}
1693  aftsTrack=anEvent->GetTrack(i2);
1694  if(!(aftsTrack->InRPSelection())){continue;}
1695  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1696  dPhi2=aftsTrack->Phi();
1697  if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1698  for(Int_t i3=0;i3<nPrim;i3++)
1699  {
1700  if(i3==i1||i3==i2){continue;}
1701  aftsTrack=anEvent->GetTrack(i3);
1702  if(!(aftsTrack->InRPSelection())){continue;}
1703  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1704  dPhi3=aftsTrack->Phi();
1705  if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1706  for(Int_t i4=0;i4<nPrim;i4++)
1707  {
1708  if(i4==i1||i4==i2||i4==i3){continue;}
1709  aftsTrack=anEvent->GetTrack(i4);
1710  if(!(aftsTrack->InRPSelection())){continue;}
1711  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1712  dPhi4=aftsTrack->Phi();
1713  if(fUseWeights[0][0]){wPhi4=Weight(dPhi4,"RP","phi");}
1714  for(Int_t i5=0;i5<nPrim;i5++)
1715  {
1716  if(i5==i1||i5==i2||i5==i3||i5==i4){continue;}
1717  aftsTrack=anEvent->GetTrack(i5);
1718  if(!(aftsTrack->InRPSelection())){continue;}
1719  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1720  dPhi5=aftsTrack->Phi();
1721  if(fUseWeights[0][0]){wPhi5=Weight(dPhi5,"RP","phi");}
1722  for(Int_t i6=0;i6<nPrim;i6++)
1723  {
1724  if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5){continue;}
1725  aftsTrack=anEvent->GetTrack(i6);
1726  if(!(aftsTrack->InRPSelection())){continue;}
1727  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1728  dPhi6=aftsTrack->Phi();
1729  if(fUseWeights[0][0]){wPhi6=Weight(dPhi6,"RP","phi");}
1730  for(Int_t i7=0;i7<nPrim;i7++)
1731  {
1732  if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6){continue;}
1733  aftsTrack=anEvent->GetTrack(i7);
1734  if(!(aftsTrack->InRPSelection())){continue;}
1735  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1736  dPhi7=aftsTrack->Phi();
1737  if(fUseWeights[0][0]){wPhi7=Weight(dPhi7,"RP","phi");}
1738  for(Int_t i8=0;i8<nPrim;i8++)
1739  {
1740  if(i8==i1||i8==i2||i8==i3||i8==i4||i8==i5||i8==i6||i8==i7){continue;}
1741  aftsTrack=anEvent->GetTrack(i8);
1742  if(!(aftsTrack->InRPSelection())){continue;}
1743  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1744  dPhi8=aftsTrack->Phi();
1745  if(fUseWeights[0][0]){wPhi8=Weight(dPhi8,"RP","phi");}
1746  // Fill:
1747  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);
1748  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);
1749  } // end of for(Int_t i8=0;i8<nPrim;i8++)
1750  } // end of for(Int_t i7=0;i7<nPrim;i7++)
1751  } // end of for(Int_t i6=0;i6<nPrim;i6++)
1752  } // end of for(Int_t i5=0;i5<nPrim;i5++)
1753  } // end of for(Int_t i4=0;i4<nPrim;i4++)
1754  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1755  } // end of for(Int_t i2=0;i2<nPrim;i2++)
1756  } // end of for(Int_t i1=0;i1<nPrim;i1++)
1757  } // end of if(nPrim>=8)
1758 
1759  // *) Printout: TBI move somewhere else
1760  printf("\n cosine:");
1761  printf("\n 1-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(2));
1762  printf("\n 1-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(1));
1763  printf("\n 2-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(4));
1764  printf("\n 2-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(3));
1765  printf("\n 3-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(6));
1766  printf("\n 3-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(5));
1767  printf("\n 4-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(8));
1768  printf("\n 4-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(7));
1769  printf("\n 5-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(10));
1770  printf("\n 5-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(9));
1771  printf("\n 6-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(12));
1772  printf("\n 6-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(11));
1773  printf("\n 7-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(14));
1774  printf("\n 7-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(13));
1775  printf("\n 8-p => Q-vector: %.12f",fNestedLoopsResultsCosPro->GetBinContent(16));
1776  printf("\n 8-p => Nested loops: %.12f",fNestedLoopsResultsCosPro->GetBinContent(15));
1777 
1778  printf("\n\n sinus:");
1779  printf("\n 1-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(2));
1780  printf("\n 1-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(1));
1781  printf("\n 2-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(4));
1782  printf("\n 2-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(3));
1783  printf("\n 3-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(6));
1784  printf("\n 3-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(5));
1785  printf("\n 4-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(8));
1786  printf("\n 4-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(7));
1787  printf("\n 5-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(10));
1788  printf("\n 5-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(9));
1789  printf("\n 6-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(12));
1790  printf("\n 6-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(11));
1791  printf("\n 7-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(14));
1792  printf("\n 7-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(13));
1793  printf("\n 8-p => Q-vector: %.12f",fNestedLoopsResultsSinPro->GetBinContent(16));
1794  printf("\n 8-p => Nested loops: %.12f",fNestedLoopsResultsSinPro->GetBinContent(15));
1795 
1796  printf("\n\n");
1797 
1798 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckWithNestedLoops(AliFlowEventSimple *anEvent)
1799 
1800 //=======================================================================================================================
1801 
1803 {
1804  // Cross-check results for differential multi-particle correlations with nested loops.
1805 
1806  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckDiffWithNestedLoops(AliFlowEventSimple *anEvent)";
1807 
1808  Int_t nPrim = anEvent->NumberOfTracks();
1809  AliFlowTrackSimple *aftsTrack = NULL;
1810  Double_t dPsi1=0.,dPhi2=0.,dPhi3=0.,dPhi4=0.;
1811  Double_t wPsi1=1.,wPhi2=1.,wPhi3=1.,wPhi4=1.;
1812 
1813  Int_t cs = fCrossCheckDiffCSCOBN[0]; // cos/sin
1814 
1815  // TBI reimplement lines below in a more civilised manner:
1816  Bool_t bCrossCheck2p = kFALSE;
1817  Bool_t bCrossCheck3p = kFALSE;
1818  Bool_t bCrossCheck4p = kFALSE;
1819 
1820  if(fCrossCheckDiffCSCOBN[1] == 2){bCrossCheck2p = kTRUE;}
1821  else if(fCrossCheckDiffCSCOBN[1] == 3){bCrossCheck3p = kTRUE;}
1822  else if(fCrossCheckDiffCSCOBN[1] == 4){bCrossCheck4p = kTRUE;}
1823 
1824  if(Int_t(bCrossCheck2p + bCrossCheck3p + bCrossCheck4p) > 1)
1825  {
1826  Fatal(sMethodName.Data(),"Int_t(bCrossCheck2p + bCrossCheck3p + bCrossCheck4p) > 1");
1827  }
1828  if(!(bCrossCheck2p || bCrossCheck3p || bCrossCheck4p))
1829  {
1830  Fatal(sMethodName.Data(),"!(bCrossCheck2p || bCrossCheck3p || bCrossCheck4p)");
1831  }
1832  Int_t nDiffBinNo = fCrossCheckDiffCSCOBN[2];
1833  Double_t dPt = 0., dEta = 0.;
1834 
1835  // <2'>:
1836  for(Int_t i1=0;i1<nPrim;i1++) // Loop over particles in a differential bin
1837  {
1838  aftsTrack=anEvent->GetTrack(i1);
1839  if(!(aftsTrack->InPOISelection())){continue;}
1840  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1841  dPsi1=aftsTrack->Phi();
1843  {
1844  dPt=aftsTrack->Pt();
1845  if(fDiffCorrelationsPro[0][1]->FindBin(dPt) != nDiffBinNo){continue;} // TBI spaghetti again
1846  } else
1847  {
1848  dEta=aftsTrack->Eta();
1849  if(fDiffCorrelationsPro[0][1]->FindBin(dEta) != nDiffBinNo){continue;} // TBI spaghetti again
1850  }
1851  if(fUseWeights[1][0]){wPsi1=Weight(dPsi1,"POI","phi");}
1852  for(Int_t i2=0;i2<nPrim;i2++) // Loop over particles in an event
1853  {
1854  if(i2==i1){continue;} // get rid of autocorrelations
1855  aftsTrack=anEvent->GetTrack(i2);
1856  if(!(aftsTrack->InRPSelection())){continue;}
1857  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1858  dPhi2=aftsTrack->Phi();
1859  if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1860  // Fill profiles:
1861  if(bCrossCheck2p)
1862  {
1863  if(fCrossCheckDiffCSCOBN[0] == 0)
1864  {
1865  fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Cos(fDiffHarmonics[1][0]*dPsi1+fDiffHarmonics[1][1]*dPhi2),wPsi1*wPhi2);
1866  } else {fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Sin(fDiffHarmonics[1][0]*dPsi1+fDiffHarmonics[1][1]*dPhi2),wPsi1*wPhi2);}
1867  } // if(bCrossCheck2p)
1868  } // for(Int_t i2=0;i2<nPrim;i2++)
1869  } // for(Int_t i1=0;i1<nPrim;i1++)
1870 
1871  // <3'>:
1872  for(Int_t i1=0;i1<nPrim;i1++) // Loop over particles in a differential bin
1873  {
1874  aftsTrack=anEvent->GetTrack(i1);
1875  if(!(aftsTrack->InPOISelection())){continue;}
1876  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1877  dPsi1=aftsTrack->Phi();
1879  {
1880  dPt=aftsTrack->Pt();
1881  if(fDiffCorrelationsPro[0][1]->FindBin(dPt) != nDiffBinNo){continue;} // TBI spaghetti again
1882  } else
1883  {
1884  dEta=aftsTrack->Eta();
1885  if(fDiffCorrelationsPro[0][1]->FindBin(dEta) != nDiffBinNo){continue;} // TBI spaghetti again
1886  }
1887  if(fUseWeights[1][0]){wPsi1=Weight(dPsi1,"POI","phi");}
1888  for(Int_t i2=0;i2<nPrim;i2++) // Loop over particles in an event
1889  {
1890  if(i2==i1){continue;} // get rid of autocorrelations
1891  aftsTrack=anEvent->GetTrack(i2);
1892  if(!(aftsTrack->InRPSelection())){continue;}
1893  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1894  dPhi2=aftsTrack->Phi();
1895  if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1896  for(Int_t i3=0;i3<nPrim;i3++)
1897  {
1898  if(i3==i1||i3==i2){continue;} // get rid of autocorrelations
1899  aftsTrack=anEvent->GetTrack(i3);
1900  if(!(aftsTrack->InRPSelection())){continue;}
1901  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1902  dPhi3=aftsTrack->Phi();
1903  if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1904  // Fill the profiles:
1905  if(bCrossCheck3p)
1906  {
1907  if(fCrossCheckDiffCSCOBN[0] == 0)
1908  {
1909  fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Cos(fDiffHarmonics[2][0]*dPsi1+fDiffHarmonics[2][1]*dPhi2+fDiffHarmonics[2][2]*dPhi3),wPsi1*wPhi2*wPhi3);
1910  } else {fNestedLoopsDiffResultsPro->Fill(0.5,TMath::Sin(fDiffHarmonics[2][0]*dPsi1+fDiffHarmonics[2][1]*dPhi2+fDiffHarmonics[2][2]*dPhi3),wPsi1*wPhi2*wPhi3);}
1911  } // if(bCrossCheck3p)
1912  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1913  } // for(Int_t i2=0;i2<nPrim;i2++)
1914  } // for(Int_t i1=0;i1<nPrim;i1++)
1915 
1916  // <4'>:
1917  for(Int_t i1=0;i1<nPrim;i1++) // Loop over particles in a differential bin
1918  {
1919  aftsTrack=anEvent->GetTrack(i1);
1920  if(!(aftsTrack->InPOISelection())){continue;}
1921  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1922  dPsi1=aftsTrack->Phi();
1924  {
1925  dPt=aftsTrack->Pt();
1926  if(fDiffCorrelationsPro[0][1]->FindBin(dPt) != nDiffBinNo){continue;} // TBI spaghetti again
1927  } else
1928  {
1929  dEta=aftsTrack->Eta();
1930  if(fDiffCorrelationsPro[0][1]->FindBin(dEta) != nDiffBinNo){continue;} // TBI spaghetti again
1931  }
1932  if(fUseWeights[1][0]){wPsi1=Weight(dPsi1,"POI","phi");}
1933  for(Int_t i2=0;i2<nPrim;i2++) // Loop over particles in an event
1934  {
1935  if(i2==i1){continue;} // get rid of autocorrelations
1936  aftsTrack=anEvent->GetTrack(i2);
1937  if(!(aftsTrack->InRPSelection())){continue;}
1938  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1939  dPhi2=aftsTrack->Phi();
1940  if(fUseWeights[0][0]){wPhi2=Weight(dPhi2,"RP","phi");}
1941  for(Int_t i3=0;i3<nPrim;i3++)
1942  {
1943  if(i3==i1||i3==i2){continue;} // get rid of autocorrelations
1944  aftsTrack=anEvent->GetTrack(i3);
1945  if(!(aftsTrack->InRPSelection())){continue;}
1946  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1947  dPhi3=aftsTrack->Phi();
1948  if(fUseWeights[0][0]){wPhi3=Weight(dPhi3,"RP","phi");}
1949  for(Int_t i4=0;i4<nPrim;i4++)
1950  {
1951  if(i4==i1||i4==i2||i4==i3){continue;} // get rid of autocorrelations
1952  aftsTrack=anEvent->GetTrack(i4);
1953  if(!(aftsTrack->InRPSelection())){continue;}
1954  if(fSkipSomeIntervals && !TrackIsInSpecifiedIntervals(aftsTrack)){continue;} // TBI tmp gym
1955  dPhi4=aftsTrack->Phi();
1956  if(fUseWeights[0][0]){wPhi4=Weight(dPhi4,"RP","phi");}
1957  // Fill the profiles:
1958  if(bCrossCheck4p)
1959  {
1960  if(fCrossCheckDiffCSCOBN[0] == 0)
1961  {
1962  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);
1963  } 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);}
1964  } // if(bCrossCheck4p)
1965  } // end of for(Int_t i4=0;i4<nPrim;i4++)
1966  } // end of for(Int_t i3=0;i3<nPrim;i3++)
1967  } // for(Int_t i2=0;i2<nPrim;i2++)
1968  } // for(Int_t i1=0;i1<nPrim;i1++)
1969 
1970  // Printout:
1971  // 2-p:
1972  if(bCrossCheck2p)
1973  {
1974  printf("\n 2-p => Q-vector: %.12f",fDiffCorrelationsPro[cs][1]->GetBinContent(nDiffBinNo));
1975  printf("\n 2-p => Nested loops: %.12f\n",fNestedLoopsDiffResultsPro->GetBinContent(1));
1976  }
1977  // 3-p:
1978  if(bCrossCheck3p)
1979  {
1980  printf("\n 3-p => Q-vector: %.12f",fDiffCorrelationsPro[cs][2]->GetBinContent(nDiffBinNo));
1981  printf("\n 3-p => Nested loops: %.12f\n",fNestedLoopsDiffResultsPro->GetBinContent(1));
1982  }
1983  // 4-p:
1984  if(bCrossCheck4p)
1985  {
1986  printf("\n 4-p => Q-vector: %.12f",fDiffCorrelationsPro[cs][3]->GetBinContent(nDiffBinNo));
1987  printf("\n 4-p => Nested loops: %.12f\n",fNestedLoopsDiffResultsPro->GetBinContent(1));
1988  }
1989 
1990 } // void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckDiffWithNestedLoops(AliFlowEventSimple *anEvent)
1991 
1992 //=======================================================================================================================
1993 
1995 {
1996  // Fill Q-vector components.
1997 
1998  Int_t nTracks = anEvent->NumberOfTracks(); // TBI shall I promote this to data member?
1999  Double_t dPhi = 0., wPhi = 1.; // azimuthal angle and corresponding phi weight
2000  Double_t dPt = 0., wPt = 1.; // transverse momentum and corresponding pT weight
2001  Double_t dEta = 0., wEta = 1.; // pseudorapidity and corresponding eta weight
2002  Double_t wToPowerP = 1.; // weight raised to power p
2003  Int_t nCounterRPs = 0;
2004  for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2005  {
2006  AliFlowTrackSimple *pTrack = NULL;
2007  if(!fSelectRandomlyRPs) // TBI hw RPs
2008  {
2009  pTrack = anEvent->GetTrack(t);
2010  }
2011  else
2012  {
2013  pTrack = anEvent->GetTrack((Int_t)fRandomIndicesRPs->GetAt(t));
2014  }
2015 
2016  if(!pTrack){printf("\n AAAARGH: pTrack is NULL in MPC::FillQvector(...) !!!!"); continue;}
2017 
2018  if(!TrackIsInSpecifiedIntervals(pTrack)){continue;} // TBI tmp gym
2019 
2020  if(!(pTrack->InRPSelection() || pTrack->InPOISelection())){printf("\n AAAARGH: pTrack is neither RP nor POI !!!!"); continue;}
2021 
2022  if(pTrack->InRPSelection()) // fill Q-vector components only with reference particles
2023  {
2024  nCounterRPs++;
2025  if(fSelectRandomlyRPs && nCounterRPs == fnSelectedRandomlyRPs){break;} // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2026 
2027  wPhi = 1.; wPt = 1.; wEta = 1.; wToPowerP = 1.; // TBI this shall go somewhere else, for performance sake
2028 
2029  // Access kinematic variables for RP and corresponding weights:
2030  dPhi = pTrack->Phi(); // azimuthal angle
2031  if(fUseWeights[0][0]){wPhi = Weight(dPhi,"RP","phi");} // corresponding phi weight
2032  //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
2033  //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
2034  dPt = pTrack->Pt();
2035  if(fUseWeights[0][1]){wPt = Weight(dPt,"RP","pt");} // corresponding pT weight
2036  dEta = pTrack->Eta();
2037  if(fUseWeights[0][2]){wEta = Weight(dEta,"RP","eta");} // corresponding eta weight
2038 
2039  // Calculate Q-vector components:
2040  for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
2041  {
2042  for(Int_t wp=0;wp<fMaxCorrelator+1;wp++) // weight power
2043  {
2044  if(fUseWeights[0][0]||fUseWeights[0][1]||fUseWeights[0][2]){wToPowerP = pow(wPhi*wPt*wEta,wp);}
2045  fQvector[h][wp] += TComplex(wToPowerP*TMath::Cos(h*dPhi),wToPowerP*TMath::Sin(h*dPhi));
2046  } // for(Int_t wp=0;wp<fMaxCorrelator+1;wp++)
2047  } // for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
2048  } // if(pTrack->InRPSelection()) // fill Q-vector components only with reference particles
2049 
2050  // Differential Q-vectors (a.k.a. p-vector and q-vector):
2051  if(!fCalculateDiffQvectors){continue;}
2052  if(pTrack->InPOISelection())
2053  {
2054  wPhi = 1.; wPt = 1.; wEta = 1.; wToPowerP = 1.; // TBI this shall go somewhere else, for performance sake
2055 
2056  // Access kinematic variables for POI and corresponding weights:
2057  dPhi = pTrack->Phi(); // azimuthal angle
2058  if(fUseWeights[1][0]){wPhi = Weight(dPhi,"POI","phi");} // corresponding phi weight
2059  //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
2060  //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
2061  dPt = pTrack->Pt();
2062  if(fUseWeights[1][1]){wPt = Weight(dPt,"POI","pt");} // corresponding pT weight
2063  dEta = pTrack->Eta();
2064  if(fUseWeights[1][2]){wEta = Weight(dEta,"POI","eta");} // corresponding eta weight
2065 
2066  // Determine bin:
2067  Int_t binNo = -44;
2069  {
2070  binNo = fDiffCorrelationsPro[0][0]->FindBin(dPt); // TBI: hardwired [0][0]
2071  } else
2072  {
2073  binNo = fDiffCorrelationsPro[0][0]->FindBin(dEta); // TBI: hardwired [0][0]
2074  }
2075  // Calculate p-vector components:
2076  for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
2077  {
2078  for(Int_t wp=0;wp<fMaxCorrelator+1;wp++) // weight power
2079  {
2080  if(fUseWeights[1][0]||fUseWeights[1][1]||fUseWeights[1][2]){wToPowerP = pow(wPhi*wPt*wEta,wp);}
2081  fpvector[binNo-1][h][wp] += TComplex(wToPowerP*TMath::Cos(h*dPhi),wToPowerP*TMath::Sin(h*dPhi));
2082 
2083  if(pTrack->InRPSelection())
2084  {
2085  // Fill q-vector components:
2086  wPhi = 1.; wPt = 1.; wEta = 1.; wToPowerP = 1.; // TBI this shall go somewhere else, for performance sake
2087 
2088  if(fUseWeights[0][0]){wPhi = Weight(dPhi,"RP","phi");} // corresponding phi weight
2089  //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
2090  //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
2091  if(fUseWeights[0][1]){wPt = Weight(dPt,"RP","pt");} // corresponding pT weight
2092  if(fUseWeights[0][2]){wEta = Weight(dEta,"RP","eta");} // corresponding eta weight
2093  if(fUseWeights[1][0]){wPhi = Weight(dPhi,"POI","phi");} // corresponding phi weight
2094  //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
2095  //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
2096  if(fUseWeights[1][1]){wPt = Weight(dPt,"POI","pt");} // corresponding pT weight
2097  if(fUseWeights[1][2]){wEta = Weight(dEta,"POI","eta");} // corresponding eta weight
2098  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);}
2099  fqvector[binNo-1][h][wp] += TComplex(wToPowerP*TMath::Cos(h*dPhi),wToPowerP*TMath::Sin(h*dPhi));
2100  } // if(pTrack->InRPSelection())
2101 
2102  } // for(Int_t wp=0;wp<fMaxCorrelator+1;wp++)
2103  } // for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
2104  } // if(pTrack->InPOISelection())
2105 
2106  } // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2107 
2108 } // void AliFlowAnalysisWithMultiparticleCorrelations::FillQvector(AliFlowEventSimple *anEvent)
2109 
2110 //=======================================================================================================================
2111 
2113 {
2114  // Cross-check all initial settings in this method.
2115 
2116  // a) Few cross-checks for control histograms;
2117  // b) Few cross-checks for flags for correlations;
2118  // c) 'Standard candles';
2119  // d) Q-cumulants;
2120  // e) Weights;
2121  // f) Differential correlations;
2122  // g) Nested loops;
2123  // h) Dump the points.
2124 
2125  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckSettings()";
2126 
2127  // a) Few cross-checks for control histograms: TBI the lines below are not really what they are supposed to be...
2128  /*
2129  if(fFillKinematicsHist && !fFillControlHistograms){Fatal(sMethodName.Data(),"fFillKinematicsHist && !fFillControlHistograms");}
2130  if(fFillMultDistributionsHist && !fFillControlHistograms){Fatal(sMethodName.Data(),"fFillMultDistributionsHist && !fFillControlHistograms");}
2131  if(fFillMultCorrelationsHist && !fFillControlHistograms){Fatal(sMethodName.Data(),"fFillMultCorrelationsHist && !fFillControlHistograms");}
2132  */
2133 
2134  // b) Few cross-checks for flags for correlations: // TBI the lines bellow can be civilized
2136  if(iSum>1){Fatal(sMethodName.Data(),"iSum is doing crazy things...");}
2137  if(fCalculateOnlyCos && fCalculateOnlySin){Fatal(sMethodName.Data(),"fCalculateOnlyCos && fCalculateOnlySin");}
2138 
2139  // c) 'Standard candles':
2141  {
2142  Fatal(sMethodName.Data(),"fCalculateStandardCandles && !fCalculateCorrelations");
2143  }
2145  {
2146  Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateCorrelations && fCalculateSameIsotropic");
2147  }
2149  {
2150  Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateOnlyForHarmonicQC");
2151  }
2153  {
2154  Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateOnlyForSC && (4!=fDontGoBeyond)");
2155  }
2157  {
2158  Warning(sMethodName.Data(),"fCalculateStandardCandles && !fPropagateErrorSC");
2159  }
2161  {
2162  Fatal(sMethodName.Data(),"fCalculateStandardCandles && fCalculateOnlySin");
2163  }
2165  {
2166  Fatal(sMethodName.Data(),"fCalculateStandardCandles && fDontGoBeyond < 3");
2167  }
2168 
2169  // d) Q-cumulants:
2171  {
2172  Fatal(sMethodName.Data(),"fCalculateQcumulants && !fCalculateCorrelations");
2173  }
2174  if(fCalculateQcumulants && !(fHarmonicQC > 0))
2175  {
2176  Fatal(sMethodName.Data(),"fCalculateQcumulants && !(fHarmonicQC > 0)");
2177  }
2179  {
2180  Fatal(sMethodName.Data(),"fCalculateQcumulants && fCalculateOnlyForSC");
2181  }
2183  {
2184  Warning(sMethodName.Data(),"fCalculateQcumulants && !fPropagateErrorQC");
2185  }
2187  {
2188  Fatal(sMethodName.Data(),"fCalculateQcumulants && fCalculateOnlySin");
2189  }
2190 
2191  // e) Weights:
2192  for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2193  {
2194  for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2195  {
2196  if(fUseWeights[rp][ppe] && !fWeightsHist[rp][ppe])
2197  {
2198  Fatal(sMethodName.Data(),"fUseWeights[rp][ppe] && !fWeightsHist[rp][ppe], rp = %d, ppe = %d",rp,ppe);
2199  }
2200  }
2201  }
2202 
2203  // f) Differential correlations:
2205  {
2206  Fatal(sMethodName.Data(),"fCalculateDiffCorrelations && !fCalculateDiffQvectors");
2207  }
2209  {
2210  Fatal(sMethodName.Data(),"fCalculateDiffCorrelations && !fCalculateQvector");
2211  }
2213  {
2214  Fatal(sMethodName.Data(),"!fCalculateDiffCorrelations && fCalculateDiffQvectors");
2215  }
2217  {
2218  Fatal(sMethodName.Data(),"fCalculateDiffCorrelations && !fUseDefaultBinning && (fnDiffBins < 1 || !fRangesDiffBins)");
2219  }
2221  {
2222  Fatal(sMethodName.Data(),"fCalculateDiffCorrelations && !(fCalculateDiffCos || fCalculateDiffSin)");
2223  }
2225  {
2226  Warning(sMethodName.Data(),"fCalculateDiffCorrelations && fDontFill[1]");
2227  }
2228 
2229  // g) Nested loops:
2231  {
2232  Fatal(sMethodName.Data(),"fCrossCheckDiffWithNestedLoops && (1 == fCrossCheckDiffCSCOBN[0] && !CalculateDiffSin)");
2233  }
2235  {
2236  Fatal(sMethodName.Data(),"fCrossCheckDiffWithNestedLoops && (0 == fCrossCheckDiffCSCOBN[0] && !CalculateDiffCos)");
2237  }
2238 
2239  // h) Dump the points:
2241  {
2242  Fatal(sMethodName.Data(),"if(fDumpThePoints && !fFillMultDistributionsHist)");
2243  }
2245  {
2246  Fatal(sMethodName.Data(),"if(fDumpThePoints && fMaxNoEventsPerFile <= 0)");
2247  }
2248 
2249 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::CrossCheckSettings()
2250 
2251 //=======================================================================================================================
2252 
2254 {
2255  // Book and nest all lists nested in the base list fHistList.
2256 
2257  // a) Book and nest lists for control histograms;
2258  // b) Book and nest lists for Q-vectors;
2259  // c) Book and nest lists for correlations;
2260  // d) Book and nest lists for e-b-e cumulants;
2261  // e) Book and nest lists for weights;
2262  // f) Book and nest lists for nested loops;
2263  // g) Book and nest lists for 'standard candles';
2264  // h) Book and nest lists for Q-cumulants;
2265  // i) Book and nest lists for differential correlations;
2266  // j) Book and nest lists for symmetry plane correlations;
2267  // k) Book and nest lists for correlations with eta gaps.
2268 
2269  // a) Book and nest lists for control histograms:
2270  fControlHistogramsList = new TList();
2271  fControlHistogramsList->SetName("Control Histograms");
2272  fControlHistogramsList->SetOwner(kTRUE);
2274 
2275  // b) Book and nest lists for Q-vectors:
2276  fQvectorList = new TList();
2277  fQvectorList->SetName("Q-vectors");
2278  fQvectorList->SetOwner(kTRUE);
2279  fHistList->Add(fQvectorList);
2280 
2281  // c) Book and nest lists for correlations:
2282  fCorrelationsList = new TList();
2283  fCorrelationsList->SetName("Correlations");
2284  fCorrelationsList->SetOwner(kTRUE);
2286 
2287  // d) Book and nest lists for e-b-e cumulants:
2288  fEbECumulantsList = new TList();
2289  fEbECumulantsList->SetName("E-b-e Cumulants");
2290  fEbECumulantsList->SetOwner(kTRUE);
2292 
2293  // e) Book and nest lists for weights:
2294  fWeightsList = new TList();
2295  fWeightsList->SetName("Weights");
2296  fWeightsList->SetOwner(kTRUE);
2297  fHistList->Add(fWeightsList);
2298 
2299  // f) Book and nest lists for nested loops:
2300  fNestedLoopsList = new TList();
2301  fNestedLoopsList->SetName("Nested Loops");
2302  fNestedLoopsList->SetOwner(kTRUE);
2304 
2305  // g) Book and nest lists for 'standard candles':
2306  fStandardCandlesList = new TList();
2307  fStandardCandlesList->SetName("Standard Candles");
2308  fStandardCandlesList->SetOwner(kTRUE);
2310 
2311  // h) Book and nest lists for Q-cumulants:
2312  fQcumulantsList = new TList();
2313  fQcumulantsList->SetName("Q-cumulants");
2314  fQcumulantsList->SetOwner(kTRUE);
2315  fHistList->Add(fQcumulantsList);
2316 
2317  // i) Book and nest lists for differential correlations:
2318  fDiffCorrelationsList = new TList();
2319  fDiffCorrelationsList->SetName("Differential Correlations");
2320  fDiffCorrelationsList->SetOwner(kTRUE);
2322 
2323  // j) Book and nest lists for symmetry plane correlations:
2324  fSymmetryPlanesList = new TList();
2325  fSymmetryPlanesList->SetName("Symmetry_Plane_Correlations");
2326  fSymmetryPlanesList->SetOwner(kTRUE);
2328 
2329  // k) Book and nest lists for correlations with eta gaps:
2330  fEtaGapsList = new TList();
2331  fEtaGapsList->SetName("Correlations_with_eta_gaps");
2332  fEtaGapsList->SetOwner(kTRUE);
2333  fHistList->Add(fEtaGapsList);
2334 
2335 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookAndNestAllLists()
2336 
2337 //=======================================================================================================================
2338 
2340 {
2341  // Store the final results in output file <outputFileName>.root.
2342 
2343  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
2344  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
2345 
2346  delete output;
2347 
2348 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::WriteHistograms(TString outputFileName)
2349 
2350 //=======================================================================================================================
2351 
2353 {
2354  // Store the final results in output file <outputFileName>.root.
2355 
2356  outputFileName->Add(fHistList);
2357  outputFileName->Write(outputFileName->GetName(),TObject::kSingleKey);
2358 
2359 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::WriteHistograms(TDirectoryFile *outputFileName)
2360 
2361 //=======================================================================================================================
2362 
2364 {
2365  // Book all the stuff for control histograms.
2366 
2367  // a) Book the profile holding all the flags for control histograms;
2368  // b) Book all control histograms;
2369  // b0) Book TH1D *fKinematicsHist[2][3];
2370  // b1) Book TH1D *fMultDistributionsHist[3];
2371  // b2) Book TH2D *fMultCorrelationsHist[3].
2372 
2373  // a) Book the profile holding all the flags for control histograms: TBI stil incomplete
2374  fControlHistogramsFlagsPro = new TProfile("fControlHistogramsFlagsPro","Flags and settings for control histograms",7,0,7);
2375  fControlHistogramsFlagsPro->SetTickLength(-0.01,"Y");
2376  fControlHistogramsFlagsPro->SetMarkerStyle(25);
2377  fControlHistogramsFlagsPro->SetLabelSize(0.04);
2378  fControlHistogramsFlagsPro->SetLabelOffset(0.02,"Y");
2379  fControlHistogramsFlagsPro->SetStats(kFALSE);
2380  fControlHistogramsFlagsPro->SetFillColor(kGray);
2381  fControlHistogramsFlagsPro->SetLineColor(kBlack);
2382  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(1,"fFillControlHistograms"); fControlHistogramsFlagsPro->Fill(0.5,fFillControlHistograms);
2383  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(2,"fFillKinematicsHist"); fControlHistogramsFlagsPro->Fill(1.5,fFillKinematicsHist);
2384  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(3,"fFillMultDistributionsHist"); fControlHistogramsFlagsPro->Fill(2.5,fFillMultDistributionsHist);
2385  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(4,"fFillMultCorrelationsHist"); fControlHistogramsFlagsPro->Fill(3.5,fFillMultCorrelationsHist);
2386  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(5,"fDontFill[0=RP]"); fControlHistogramsFlagsPro->Fill(4.5,fDontFill[0]);
2387  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(6,"fDontFill[1=POI]"); fControlHistogramsFlagsPro->Fill(5.5,fDontFill[1]);
2388  fControlHistogramsFlagsPro->GetXaxis()->SetBinLabel(7,"fDontFill[2=REF]"); fControlHistogramsFlagsPro->Fill(6.5,fDontFill[2]);
2390 
2391  if(!fFillControlHistograms){return;} // TBI is this safe? Well, perhaps it is if I can't implement it better...
2392 
2393  // b) Book all control histograms: // TBI add setters for all these values
2394  // b0) Book TH1D *fKinematicsHist[2][3]:
2395  TString name[2][3] = {{"RP,phi","RP,pt","RP,eta"},{"POI,phi","POI,pt","POI,eta"}}; // [RP,POI][phi,pt,eta]
2396  TString title[2] = {"Reference particles (RP)","Particles of interest (POI)"}; // [RP,POI]
2397  Int_t lineColor[2] = {kBlue,kRed}; // [RP,POI]
2398  Int_t fillColor[2] = {kBlue-10,kRed-10}; // [RP,POI]
2399  TString xAxisTitle[3] = {"#phi","p_{T}","#eta"}; // [phi,pt,eta]
2401  {
2402  for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2403  {
2404  if(fDontFill[rp]){continue;}
2405  for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2406  {
2407  fKinematicsHist[rp][ppe] = new TH1D(name[rp][ppe].Data(),title[rp].Data(),fnBins[rp][ppe],fMin[rp][ppe],fMax[rp][ppe]);
2408  fKinematicsHist[rp][ppe]->GetXaxis()->SetTitle(xAxisTitle[ppe].Data());
2409  fKinematicsHist[rp][ppe]->SetLineColor(lineColor[rp]);
2410  fKinematicsHist[rp][ppe]->SetFillColor(fillColor[rp]);
2411  fKinematicsHist[rp][ppe]->SetMinimum(0.);
2412  fControlHistogramsList->Add(fKinematicsHist[rp][ppe]);
2413  }
2414  }
2415  } // if(fFillKinematicsHist)
2416 
2417  // b1) Book TH1D *fMultDistributionsHist[3]: // TBI add setters for all these values
2418  TString nameMult[3] = {"Multiplicity (RP)","Multiplicity (POI)","Multiplicity (REF)"}; // [RP,POI,reference multiplicity]
2419  TString titleMult[3] = {"Reference particles (RP)","Particles of interest (POI)",""}; // [RP,POI,reference multiplicity]
2420  Int_t lineColorMult[3] = {kBlue,kRed,kGreen+2}; // [RP,POI,reference multiplicity]
2421  Int_t fillColorMult[3] = {kBlue-10,kRed-10,kGreen-10}; // [RP,POI,reference multiplicity]
2422  TString xAxisTitleMult[3] = {"Multiplicity (RP)","Multiplicity (POI)","Multiplicity (REF)"}; // [phi,pt,eta]
2424  {
2425  for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2426  {
2427  if(fDontFill[rprm]){continue;}
2428  fMultDistributionsHist[rprm] = new TH1D(nameMult[rprm].Data(),titleMult[rprm].Data(),fnBinsMult[rprm],fMinMult[rprm],fMaxMult[rprm]);
2429  fMultDistributionsHist[rprm]->GetXaxis()->SetTitle(xAxisTitleMult[rprm].Data());
2430  fMultDistributionsHist[rprm]->SetLineColor(lineColorMult[rprm]);
2431  fMultDistributionsHist[rprm]->SetFillColor(fillColorMult[rprm]);
2433  } // for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2434  } // if(fFillMultDistributionsHist)
2435 
2436  // b2) Book TH2I *fMultCorrelationsHist[3]:
2438  {
2439  if(!fDontFill[0] && !fDontFill[1])
2440  {
2441  fMultCorrelationsHist[0] = new TH2I("Multiplicity (RP vs. POI)","Multiplicity (RP vs. POI)",fnBinsMult[0],fMinMult[0],fMaxMult[0],fnBinsMult[1],fMinMult[1],fMaxMult[1]);
2442  fMultCorrelationsHist[0]->GetXaxis()->SetTitle(xAxisTitleMult[0].Data());
2443  fMultCorrelationsHist[0]->GetYaxis()->SetTitle(xAxisTitleMult[1].Data());
2445  }
2446  if(!fDontFill[0] && !fDontFill[2])
2447  {
2448  fMultCorrelationsHist[1] = new TH2I("Multiplicity (RP vs. REF)","Multiplicity (RP vs. REF)",fnBinsMult[0],fMinMult[0],fMaxMult[0],fnBinsMult[2],fMinMult[2],fMaxMult[2]);
2449  fMultCorrelationsHist[1]->GetXaxis()->SetTitle(xAxisTitleMult[0].Data());
2450  fMultCorrelationsHist[1]->GetYaxis()->SetTitle(xAxisTitleMult[2].Data());
2452  }
2453  if(!fDontFill[1] && !fDontFill[2])
2454  {
2455  fMultCorrelationsHist[2] = new TH2I("Multiplicity (POI vs. REF)","Multiplicity (POI vs. REF)",fnBinsMult[1],fMinMult[1],fMaxMult[1],fnBinsMult[2],fMinMult[2],fMaxMult[2]);
2456  fMultCorrelationsHist[2]->GetXaxis()->SetTitle(xAxisTitleMult[1].Data());
2457  fMultCorrelationsHist[2]->GetYaxis()->SetTitle(xAxisTitleMult[2].Data());
2459  }
2460  } // if(fFillMultCorrelationsHist){
2461 
2462 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForControlHistograms()
2463 
2464 //=======================================================================================================================
2465 
2467 {
2468  // Fill control histograms.
2469  // a) Fill TH1D *fKinematicsHist[2][3];
2470  // b) Fill TH1D *fMultDistributionsHist[3];
2471  // c) Fill TH2D *fMultCorrelationsHist[3].
2472 
2473  // a) Fill TH1D *fKinematicsHist[2][3]:
2474  Int_t nCounterRPs = 0;
2476  {
2477  Int_t nTracks = anEvent->NumberOfTracks(); // TBI shall I promote this to data member?
2478  for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2479  {
2480  AliFlowTrackSimple *pTrack = NULL;
2481  if(!fSelectRandomlyRPs) // TBI hw RPs
2482  {
2483  pTrack = anEvent->GetTrack(t);
2484  }
2485  else
2486  {
2487  pTrack = anEvent->GetTrack((Int_t)fRandomIndicesRPs->GetAt(t));
2488  }
2489  if(!pTrack){printf("\n AAAARGH: pTrack is NULL in MPC::FCH() !!!!");continue;}
2490  if(pTrack)
2491  {
2492 
2493  if(!TrackIsInSpecifiedIntervals(pTrack)){continue;} // TBI tmp gym
2494 
2495  if(pTrack->InRPSelection())
2496  {
2497  nCounterRPs++;
2498  if(fSelectRandomlyRPs && nCounterRPs == fnSelectedRandomlyRPs){break;} // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2499  }
2500 
2501  Double_t dPhi = pTrack->Phi();
2502  //if(dPhi > TMath::TwoPi()){dPhi -= TMath::TwoPi();} TBI
2503  //if(dPhi < 0.){dPhi += TMath::TwoPi();} TBI
2504  Double_t dPt = pTrack->Pt();
2505  Double_t dEta = pTrack->Eta();
2506  Double_t dPhiPtEta[3] = {dPhi,dPt,dEta};
2507 
2508  for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2509  {
2510  for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2511  {
2512  if((0==rp && pTrack->InRPSelection()) || (1==rp && pTrack->InPOISelection())) // TBI
2513  {
2514  if(fKinematicsHist[rp][ppe]){fKinematicsHist[rp][ppe]->Fill(dPhiPtEta[ppe]);}
2515  }
2516  } // for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2517  } // for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2518  } // if(pTrack)
2519  } // for(Int_t t=0;t<nTracks;t++) // loop over all tracks
2520  } // if(fFillKinematicsHist)
2521 
2522  // b) Fill TH1D *fMultDistributionsHist[3]:
2523  Double_t dMultRP = fSelectRandomlyRPs ? nCounterRPs : anEvent->GetNumberOfRPs();
2524  Double_t dMultPOI = anEvent->GetNumberOfPOIs(); // TBI reimplement when reshuffling is enabled, add support also for the POIs
2525  Double_t dMultREF = anEvent->GetReferenceMultiplicity();
2526 
2527  if(fSkipSomeIntervals) // TBI tmp gym
2528  {
2529  dMultRP = dMultRP - fNumberOfSkippedRPParticles;
2530  dMultPOI = -44;
2531  }
2532 
2533  Double_t dMult[3] = {dMultRP,dMultPOI,dMultREF};
2534  for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2535  {
2536  if(fFillMultDistributionsHist && fMultDistributionsHist[rprm]){fMultDistributionsHist[rprm]->Fill(dMult[rprm]);}
2537  }
2538 
2539  // c) Fill TH2I *fMultCorrelationsHist[3]:
2541  {
2542  if(fMultCorrelationsHist[0]){fMultCorrelationsHist[0]->Fill((Int_t)dMultRP,(Int_t)dMultPOI);} // RP vs. POI
2543  if(fMultCorrelationsHist[1]){fMultCorrelationsHist[1]->Fill((Int_t)dMultRP,(Int_t)dMultREF);} // RP vs. refMult
2544  if(fMultCorrelationsHist[2]){fMultCorrelationsHist[2]->Fill((Int_t)dMultPOI,(Int_t)dMultREF);} // POI vs. refMult
2545  } // if(fFillMultCorrelationsHist)
2546 
2547 } // void AliFlowAnalysisWithMultiparticleCorrelations::FillControlHistograms(AliFlowEventSimple *anEvent)
2548 
2549 //=======================================================================================================================
2550 
2552 {
2553  // Initialize all arrays for control histograms.
2554 
2555  // a) Initialize TH1D *fKinematicsHist[2][3];
2556  // b) Initialize TH1D *fMultDistributionsHist[3];
2557  // c) Initialize TH2D *fMultCorrelationsHist[3];
2558  // d) Initialize Bool_t fDontFill[3];
2559  // e) Initialize default binning values for fKinematicsHist[2][3];
2560  // f) Initialize default binning values for fMultCorrelationsHist[3];
2561  // g) Initialize default rp, phi and eta intervals to be skipped.
2562 
2563  // a) Initialize TH1D *fKinematicsHist[2][3]:
2564  for(Int_t rp=0;rp<2;rp++) // [RP,POI]
2565  {
2566  for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2567  {
2568  fKinematicsHist[rp][ppe] = NULL;
2569  }
2570  }
2571 
2572  // b) Initialize TH1D *fMultDistributionsHist[3]:
2573  for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
2574  {
2575  fMultDistributionsHist[rprm] = NULL;
2576  }
2577 
2578  // c) Initialize TH2I *fMultCorrelationsHist[3]:
2579  for(Int_t r=0;r<3;r++) // [RP vs. POI, RP vs. refMult, POI vs. refMult]
2580  {
2581  fMultCorrelationsHist[r] = NULL;
2582  }
2583 
2584  // d) Initialize Bool_t fDontFill[3]:
2585  for(Int_t rpr=0;rpr<3;rpr++) // [RP,POI,REF]
2586  {
2587  fDontFill[rpr] = kFALSE;
2588  }
2589 
2590  // e) Initialize default binning values for fKinematicsHist[2][3]:
2591  // nBins:
2592  fnBins[0][0] = 360; // [RP][phi]
2593  fnBins[0][1] = 1000; // [RP][pt]
2594  fnBins[0][2] = 1000; // [RP][eta]
2595  fnBins[1][0] = 360; // [POI][phi]
2596  fnBins[1][1] = 1000; // [POI][pt]
2597  fnBins[1][2] = 1000; // [POI][eta]
2598  // Min:
2599  fMin[0][0] = 0.; // [RP][phi]
2600  fMin[0][1] = 0.; // [RP][pt]
2601  fMin[0][2] = -1.; // [RP][eta]
2602  fMin[1][0] = 0.; // [POI][phi]
2603  fMin[1][1] = 0.; // [POI][pt]
2604  fMin[1][2] = -1.; // [POI][eta]
2605  // Max:
2606  fMax[0][0] = TMath::TwoPi(); // [RP][phi]
2607  fMax[0][1] = 10.; // [RP][pt]
2608  fMax[0][2] = 1.; // [RP][eta]
2609  fMax[1][0] = TMath::TwoPi(); // [POI][phi]
2610  fMax[1][1] = 10.; // [POI][pt]
2611  fMax[1][2] = 1.; // [POI][eta]
2612 
2613  // f) Initialize default binning values for fMultCorrelationsHist[3]:
2614  // nBins:
2615  fnBinsMult[0] = 3000; // [RP]
2616  fnBinsMult[1] = 3000; // [POI]
2617  fnBinsMult[2] = 3000; // [REF]
2618  // Min:
2619  fMinMult[0] = 0.; // [RP]
2620  fMinMult[1] = 0.; // [POI]
2621  fMinMult[2] = 0.; // [REF]
2622  // Max:
2623  fMaxMult[0] = 3000.; // [RP]
2624  fMaxMult[1] = 3000.; // [POI]
2625  fMaxMult[2] = 3000.; // [REF]
2626 
2627  // g) Initialize default rp, phi and eta intervals to be skipped:
2628  for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
2629  {
2630  for(Int_t i=0;i<10;i++) // interval boundaries, 10 boundaries at max
2631  {
2632  fSkip[ppe][i] = -44.;
2633  }
2634  } // for(Int_t ppe=0;ppe<3;ppe++)
2635 
2636 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForControlHistograms()
2637 
2638 //=======================================================================================================================
2639 
2641 {
2642  // Book all the stuff for Q-vector.
2643 
2644  // a) Book the profile holding all the flags for Q-vector;
2645  // ...
2646 
2647  // a) Book the profile holding all the flags for Q-vector:
2648  fQvectorFlagsPro = new TProfile("fQvectorFlagsPro","Flags for Q-vectors",2,0,2);
2649  fQvectorFlagsPro->SetTickLength(-0.01,"Y");
2650  fQvectorFlagsPro->SetMarkerStyle(25);
2651  fQvectorFlagsPro->SetLabelSize(0.03);
2652  fQvectorFlagsPro->SetLabelOffset(0.02,"Y");
2653  fQvectorFlagsPro->SetStats(kFALSE);
2654  fQvectorFlagsPro->SetFillColor(kGray);
2655  fQvectorFlagsPro->SetLineColor(kBlack);
2656  fQvectorFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateQvector"); fQvectorFlagsPro->Fill(0.5,fCalculateQvector);
2657  fQvectorFlagsPro->GetXaxis()->SetBinLabel(2,"fCalculateDiffQvectors"); fQvectorFlagsPro->Fill(1.5,fCalculateDiffQvectors);
2659 
2660  // ...
2661 
2662 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForQvector()
2663 
2664 //=======================================================================================================================
2665 
2667 {
2668  // Book all the stuff for correlations.
2669 
2670  // TBI this method can be implemented in a much more civilised way.
2671 
2672  // a) Book the profile holding all the flags for correlations;
2673  // b) Book TProfile *fCorrelationsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]).
2674 
2675  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForCorrelations()";
2676 
2677  // a) Book the profile holding all the flags for correlations:
2678  fCorrelationsFlagsPro = new TProfile("fCorrelationsFlagsPro","Flags for correlations",13,0,13);
2679  fCorrelationsFlagsPro->SetTickLength(-0.01,"Y");
2680  fCorrelationsFlagsPro->SetMarkerStyle(25);
2681  fCorrelationsFlagsPro->SetLabelSize(0.03);
2682  fCorrelationsFlagsPro->SetLabelOffset(0.02,"Y");
2683  fCorrelationsFlagsPro->SetStats(kFALSE);
2684  fCorrelationsFlagsPro->SetFillColor(kGray);
2685  fCorrelationsFlagsPro->SetLineColor(kBlack);
2686  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateCorrelations"); fCorrelationsFlagsPro->Fill(0.5,fCalculateCorrelations);
2687  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(2,"fMaxHarmonic"); fCorrelationsFlagsPro->Fill(1.5,fMaxHarmonic);
2688  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(3,"fMaxCorrelator"); fCorrelationsFlagsPro->Fill(2.5,fMaxCorrelator);
2689  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(4,"fCalculateIsotropic"); fCorrelationsFlagsPro->Fill(3.5,fCalculateIsotropic);
2690  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(5,"fCalculateSame"); fCorrelationsFlagsPro->Fill(4.5,fCalculateSame);
2691  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(6,"fSkipZeroHarmonics"); fCorrelationsFlagsPro->Fill(5.5,fSkipZeroHarmonics);
2692  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(7,"fCalculateSameIsotropic"); fCorrelationsFlagsPro->Fill(6.5,fCalculateSameIsotropic);
2693  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(8,"fCalculateAll"); fCorrelationsFlagsPro->Fill(7.5,fCalculateAll);
2694  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(9,"fDontGoBeyond"); fCorrelationsFlagsPro->Fill(8.5,fDontGoBeyond);
2695  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(10,"fCalculateOnlyForHarmonicQC"); fCorrelationsFlagsPro->Fill(9.5,fCalculateOnlyForHarmonicQC);
2696  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(11,"fCalculateOnlyForSC"); fCorrelationsFlagsPro->Fill(10.5,fCalculateOnlyForSC);
2697  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(12,"fCalculateOnlyCos"); fCorrelationsFlagsPro->Fill(11.5,fCalculateOnlyCos);
2698  fCorrelationsFlagsPro->GetXaxis()->SetBinLabel(13,"fCalculateOnlySin"); fCorrelationsFlagsPro->Fill(12.5,fCalculateOnlySin);
2700 
2701  if(!fCalculateCorrelations){return;} // TBI is this safe enough?
2702 
2703  // b) Book TProfile *fCorrelationsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]): // TBI hardwired 8, shall be fMaxCorrelator
2704  cout<<" => Booking TProfile *fCorrelationsPro[2][8]..."<<endl;
2705  TString sCosSin[2] = {"Cos","Sin"};
2706  Int_t markerColor[2] = {kBlue,kRed};
2707  Int_t markerStyle[2] = {kFullSquare,kFullSquare};
2708  Int_t nBins[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
2709  Int_t nBinsTitle[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
2710  Int_t nToBeFilled[8] = {0,0,0,0,0,0,0,0}; // TBI hardwired 8, shall be fMaxCorrelator
2711  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2712  {
2713  // Implementing \binom{n+k-1}{k}, which is the resulting number of sets obtained
2714  // after sampling n starting elements into k subsets, repetitions allowed.
2715  // In my case, n=2*fMaxHarmonic+1, k=c+1, hence:
2716  nBins[c] = (Int_t)(TMath::Factorial(2*fMaxHarmonic+1+c+1-1)
2717  / (TMath::Factorial(2*fMaxHarmonic+1-1)*TMath::Factorial(c+1)));
2718  nBinsTitle[c] = nBins[c];
2719  if(c>=fDontGoBeyond){nBins[c]=1;} // TBI is this really safe?
2720  } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
2721  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2722  {
2723  if(fCalculateOnlyCos && 1==cs){continue;}
2724  else if(fCalculateOnlySin && 0==cs){continue;}
2725  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2726  {
2727  if(c==fDontGoBeyond){continue;}
2728  if(fCalculateOnlyForHarmonicQC && c%2==0){continue;}
2729  fCorrelationsPro[cs][c] = new TProfile(Form("%dpCorrelations%s",c+1,sCosSin[cs].Data()),"",nBins[c],0.,1.*nBins[c]);
2730  fCorrelationsPro[cs][c]->Sumw2();
2731  fCorrelationsPro[cs][c]->SetStats(kFALSE);
2732  fCorrelationsPro[cs][c]->SetMarkerColor(markerColor[cs]);
2733  fCorrelationsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
2734  fCorrelationsList->Add(fCorrelationsPro[cs][c]);
2735  } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
2736  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2737  // Set all bin labels: TBI this can be implemented better, most likely...
2738  Int_t binNo[2][8];
2739  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2740  {
2741  if(fCalculateOnlyCos && 1==cs){continue;}
2742  else if(fCalculateOnlySin && 0==cs){continue;}
2743  for(Int_t c=0;c<fMaxCorrelator;c++)
2744  {
2745  binNo[cs][c] = 1;
2746  }
2747  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2748 
2749  for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++)
2750  {
2751  cout<< Form(" Patience, this takes some time... n1 = %d/%d\r",n1+fMaxHarmonic,2*fMaxHarmonic)<<flush; // TBI
2752  if(fSkipZeroHarmonics && 0==n1){continue;}
2753  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n1) != fHarmonicQC){continue;}
2754  if(fCalculateAll)
2755  {
2756  for(Int_t cs=0;cs<2;cs++)
2757  {
2758  if(fCalculateOnlyCos && 1==cs){continue;}
2759  else if(fCalculateOnlySin && 0==cs){continue;}
2760  if(fCorrelationsPro[cs][0]){fCorrelationsPro[cs][0]->GetXaxis()->SetBinLabel(binNo[cs][0]++,Form("%s(%d)",sCosSin[cs].Data(),n1));}
2761  } // for(Int_t cs=0;cs<2;cs++)
2762  nToBeFilled[0]++;
2763  }
2764  if(1==fDontGoBeyond){continue;}
2765  for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
2766  {
2767  if(fSkipZeroHarmonics && 0==n2){continue;}
2768  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n2) != fHarmonicQC){continue;}
2769  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2))
2770  || (fCalculateSameIsotropic && 0==n1+n2 && TMath::Abs(n1)==TMath::Abs(n2))
2771  || (fCalculateOnlyForHarmonicQC && 0==n1+n2)
2772  || (fCalculateOnlyForSC && 0==n1+n2))
2773  {
2774  for(Int_t cs=0;cs<2;cs++)
2775  {
2776  if(fCalculateOnlyCos && 1==cs){continue;}
2777  else if(fCalculateOnlySin && 0==cs){continue;}
2778  if(fCorrelationsPro[cs][1]){fCorrelationsPro[cs][1]->GetXaxis()->SetBinLabel(binNo[cs][1]++,Form("%s(%d,%d)",sCosSin[cs].Data(),n1,n2));}
2779  } // for(Int_t cs=0;cs<2;cs++)
2780  nToBeFilled[1]++;
2781  }
2782  if(2==fDontGoBeyond){continue;}
2783  for(Int_t n3=n2;n3<=fMaxHarmonic;n3++)
2784  {
2785  if(fSkipZeroHarmonics && 0==n3){continue;}
2786  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n3) != fHarmonicQC){continue;}
2787  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3))
2788  || (fCalculateSameIsotropic && 0==n1+n2+n3 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3))
2789  || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3))
2790  {
2791  for(Int_t cs=0;cs<2;cs++)
2792  {
2793  if(fCalculateOnlyCos && 1==cs){continue;}
2794  else if(fCalculateOnlySin && 0==cs){continue;}
2795  if(fCorrelationsPro[cs][2]){fCorrelationsPro[cs][2]->GetXaxis()->SetBinLabel(binNo[cs][2]++,Form("%s(%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3));}
2796  } // for(Int_t cs=0;cs<2;cs++)
2797  nToBeFilled[2]++;
2798  }
2799  if(3==fDontGoBeyond){continue;}
2800  for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)
2801  {
2802  if(fSkipZeroHarmonics && 0==n4){continue;}
2803  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n4) != fHarmonicQC){continue;}
2804  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))
2805  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4))
2806  || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4)
2807  || (fCalculateOnlyForSC && (0==n1+n4 && 0==n2+n3 && n1 != n2 && n3 != n4)))
2808  {
2809  for(Int_t cs=0;cs<2;cs++)
2810  {
2811  if(fCalculateOnlyCos && 1==cs){continue;}
2812  else if(fCalculateOnlySin && 0==cs){continue;}
2813  if(fCorrelationsPro[cs][3]){fCorrelationsPro[cs][3]->GetXaxis()->SetBinLabel(binNo[cs][3]++,Form("%s(%d,%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3,n4));}
2814  } // for(Int_t cs=0;cs<2;cs++)
2815  nToBeFilled[3]++;
2816  }
2817  if(4==fDontGoBeyond){continue;}
2818  for(Int_t n5=n4;n5<=fMaxHarmonic;n5++)
2819  {
2820  if(fSkipZeroHarmonics && 0==n5){continue;}
2821  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n5) != fHarmonicQC){continue;}
2822  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5)
2823  || (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))
2824  || (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))
2825  || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5))
2826  {
2827  for(Int_t cs=0;cs<2;cs++)
2828  {
2829  if(fCalculateOnlyCos && 1==cs){continue;}
2830  else if(fCalculateOnlySin && 0==cs){continue;}
2831  if(fCorrelationsPro[cs][4]){fCorrelationsPro[cs][4]->GetXaxis()->SetBinLabel(binNo[cs][4]++,Form("%s(%d,%d,%d,%d,%d)",sCosSin[cs].Data(),n1,n2,n3,n4,n5));}
2832  } // for(Int_t cs=0;cs<2;cs++)
2833  nToBeFilled[4]++;
2834  }
2835  if(5==fDontGoBeyond){continue;}
2836  for(Int_t n6=n5;n6<=fMaxHarmonic;n6++)
2837  {
2838  if(fSkipZeroHarmonics && 0==n6){continue;}
2839  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n6) != fHarmonicQC){continue;}
2840  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6)
2841  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
2842  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6))
2843  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
2844  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6))
2845  || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5+n6))
2846  {
2847  for(Int_t cs=0;cs<2;cs++)
2848  {
2849  if(fCalculateOnlyCos && 1==cs){continue;}
2850  else if(fCalculateOnlySin && 0==cs){continue;}
2851  if(fCorrelationsPro[cs][5]){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));}
2852  } // for(Int_t cs=0;cs<2;cs++)
2853  nToBeFilled[5]++;
2854  }
2855  if(6==fDontGoBeyond){continue;}
2856  for(Int_t n7=n6;n7<=fMaxHarmonic;n7++)
2857  {
2858  if(fSkipZeroHarmonics && 0==n7){continue;}
2859  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n7) != fHarmonicQC){continue;}
2860  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7)
2861  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
2862  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7))
2863  || (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)
2864  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7))
2865  || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5+n6+n7))
2866  {
2867  for(Int_t cs=0;cs<2;cs++)
2868  {
2869  if(fCalculateOnlyCos && 1==cs){continue;}
2870  else if(fCalculateOnlySin && 0==cs){continue;}
2871  if(fCorrelationsPro[cs][6]){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));}
2872  } // for(Int_t cs=0;cs<2;cs++)
2873  nToBeFilled[6]++;
2874  }
2875  if(7==fDontGoBeyond){continue;}
2876  for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
2877  {
2878  if(fSkipZeroHarmonics && 0==n8){continue;}
2879  if(fCalculateOnlyForHarmonicQC && TMath::Abs(n8) != fHarmonicQC){continue;}
2880  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8)
2881  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
2882  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7) && TMath::Abs(n1)==TMath::Abs(n8))
2883  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
2884  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7)
2885  && TMath::Abs(n1)==TMath::Abs(n8))
2886  || (fCalculateOnlyForHarmonicQC && 0==n1+n2+n3+n4+n5+n6+n7+n8))
2887  {
2888  for(Int_t cs=0;cs<2;cs++)
2889  {
2890  if(fCalculateOnlyCos && 1==cs){continue;}
2891  else if(fCalculateOnlySin && 0==cs){continue;}
2892  if(fCorrelationsPro[cs][7]){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));}
2893  } // for(Int_t cs=0;cs<2;cs++)
2894  nToBeFilled[7]++;
2895  }
2896  } // for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
2897  } // for(Int_t n7=n6;n7<=fMaxHarmonic;n7++)
2898  } // for(Int_t n6=n5;n6<=fMaxHarmonic;n6++)
2899  } // for(Int_t n5=n4;n5<=fMaxHarmonic;n5++)
2900  } // for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)
2901  } // for(Int_t n3=n2;n3<=fMaxHarmonic;n3++)
2902  } // for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
2903  } // for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++)
2904 
2905  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2906  {
2907  if(fCalculateOnlyCos && 1==cs){continue;}
2908  else if(fCalculateOnlySin && 0==cs){continue;}
2909  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
2910  {
2911  if(!fCorrelationsPro[cs][c]){continue;}
2912  fCorrelationsPro[cs][c]->SetTitle(Form("%d-p correlations, %s terms, %d/%d in total",c+1,sCosSin[cs].Data(),nToBeFilled[c],nBinsTitle[c]));
2913  fCorrelationsPro[cs][c]->GetXaxis()->SetRangeUser(0.,fCorrelationsPro[cs][c]->GetBinLowEdge(nToBeFilled[c]+1));
2914  }
2915  }
2916  cout<<" Booked. "<<endl; // TBI
2917 
2918 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForCorrelations()
2919 
2920 //=======================================================================================================================
2921 
2923 {
2924  // Book all the stuff for differential correlations.
2925 
2926  // a) Book the profile holding all the flags for differential correlations;
2927  // b) Book TProfile *fDiffCorrelationsPro[2][4] ([0=cos,1=sin][1p,2p,3p,4p]).
2928 
2929  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForDiffCorrelations()";
2930 
2931  // a) Book the profile holding all the flags for differential correlations:
2932  fDiffCorrelationsFlagsPro = new TProfile("fDiffCorrelationsFlagsPro","Flags for differential correlations",5,0,5);
2933  fDiffCorrelationsFlagsPro->SetTickLength(-0.01,"Y");
2934  fDiffCorrelationsFlagsPro->SetMarkerStyle(25);
2935  fDiffCorrelationsFlagsPro->SetLabelSize(0.03);
2936  fDiffCorrelationsFlagsPro->SetLabelOffset(0.02,"Y");
2937  fDiffCorrelationsFlagsPro->SetStats(kFALSE);
2938  fDiffCorrelationsFlagsPro->SetFillColor(kGray);
2939  fDiffCorrelationsFlagsPro->SetLineColor(kBlack);
2940  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateDiffCorrelations"); fDiffCorrelationsFlagsPro->Fill(0.5,fCalculateDiffCorrelations);
2941  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(2,"fCalculateDiffCos"); fDiffCorrelationsFlagsPro->Fill(1.5,fCalculateDiffCos);
2942  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(3,"fCalculateDiffSin"); fDiffCorrelationsFlagsPro->Fill(2.5,fCalculateDiffSin);
2943  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(4,"fCalculateDiffCorrelationsVsPt"); fDiffCorrelationsFlagsPro->Fill(3.5,fCalculateDiffCorrelationsVsPt);
2944  fDiffCorrelationsFlagsPro->GetXaxis()->SetBinLabel(5,"fUseDefaultBinning"); fDiffCorrelationsFlagsPro->Fill(4.5,fUseDefaultBinning);
2946 
2947  if(!fCalculateDiffCorrelations){return;}
2948 
2949  // b) Book TProfile *fDiffCorrelationsPro[2][4] ([0=cos,1=sin][1p,2p,3p,4p]):
2950  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
2951  Int_t markerColor[2] = {kRed,kGreen};
2952  Int_t markerStyle[2] = {kFullSquare,kOpenSquare};
2953  TString sCosSin[2] = {"Cos","Sin"};
2954  TString sLabel[4] = {Form("%d",fDiffHarmonics[0][0]),
2955  Form("%d,%d",fDiffHarmonics[1][0],fDiffHarmonics[1][1]),
2956  Form("%d,%d,%d",fDiffHarmonics[2][0],fDiffHarmonics[2][1],fDiffHarmonics[2][2]),
2957  Form("%d,%d,%d,%d",fDiffHarmonics[3][0],fDiffHarmonics[3][1],fDiffHarmonics[3][2],fDiffHarmonics[3][3])};
2958 
2959  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
2960  {
2961  if(!fCalculateDiffCos && 0==cs){continue;}
2962  if(!fCalculateDiffSin && 1==cs){continue;}
2963 
2964  for(Int_t c=0;c<4;c++) // [1p,2p,3p,4p]
2965  {
2967  {
2968  if(fUseDefaultBinning)
2969  {
2970  // vs pt, default binning:
2971  fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"pt"),
2972  Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
2973  100,0.,10.);
2974  } else // if(fUseDefaultBinning)
2975  {
2976  // vs pt, non-default binning:
2977  fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"pt"),
2978  Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
2980  }// else // if(fUseDefaultBinning)
2981  fDiffCorrelationsPro[cs][c]->Sumw2();
2982  fDiffCorrelationsPro[cs][c]->SetStats(kFALSE);
2983  fDiffCorrelationsPro[cs][c]->SetMarkerColor(markerColor[cs]);
2984  fDiffCorrelationsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
2985  fDiffCorrelationsPro[cs][c]->GetXaxis()->SetTitle("p_{T}");
2986  if(fDiffStore[cs][c]){fDiffCorrelationsList->Add(fDiffCorrelationsPro[cs][c]);}
2987  } else // if(fCalculateDiffCorrelationsVsPt)
2988  {
2989  if(fUseDefaultBinning)
2990  {
2991  // vs eta, default binning:
2992  fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"eta"),
2993  Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
2994  100,-1.,1.);
2995  } else // if(fUseDefaultBinning)
2996  {
2997  // vs eta, non-default binning:
2998  fDiffCorrelationsPro[cs][c] = new TProfile(Form("%s, %dp, %s",sCosSin[cs].Data(),c+1,"eta"),
2999  Form("%s(%s)",sCosSin[cs].Data(),sLabel[c].Data()),
3001  } // else // if(fUseDefaultBinning)
3002  fDiffCorrelationsPro[cs][c]->Sumw2();
3003  fDiffCorrelationsPro[cs][c]->SetStats(kFALSE);
3004  fDiffCorrelationsPro[cs][c]->SetMarkerColor(markerColor[cs]);
3005  fDiffCorrelationsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
3006  fDiffCorrelationsPro[cs][c]->GetXaxis()->SetTitle("#eta");
3007  if(fDiffStore[cs][c]){fDiffCorrelationsList->Add(fDiffCorrelationsPro[cs][c]);}
3008  } // else // if(fCalculateDiffCorrelationsVsPt)
3009  } // for(Int_t c=0;c<4;c++) // [1p,2p,3p,4p]
3010  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
3011 
3012 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForDiffCorrelations()
3013 
3014 //=======================================================================================================================
3015 
3017 {
3018  // Book all the stuff for symmetry plane correlations.
3019 
3020  // a) Book the profile holding all the flags for symmetry plane correlations;
3021  // b) Book TProfile *fSymmetryPlanesPro[?][?]. TBI check the exact dimensions in the header file.
3022 
3023  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForSymmetryPlanes()";
3024 
3025  // a) Book the profile holding all the flags for symmetry plane correlations:
3026  fSymmetryPlanesFlagsPro = new TProfile("fSymmetryPlanesFlagsPro","Flags for Symmetry Plane Correlations (SPC)",1,0,1);
3027  fSymmetryPlanesFlagsPro->SetTickLength(-0.01,"Y");
3028  fSymmetryPlanesFlagsPro->SetMarkerStyle(25);
3029  fSymmetryPlanesFlagsPro->SetLabelSize(0.03);
3030  fSymmetryPlanesFlagsPro->SetLabelOffset(0.02,"Y");
3031  fSymmetryPlanesFlagsPro->SetStats(kFALSE);
3032  fSymmetryPlanesFlagsPro->SetFillColor(kGray);
3033  fSymmetryPlanesFlagsPro->SetLineColor(kBlack);
3034  fSymmetryPlanesFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateSymmetryPlanes"); fSymmetryPlanesFlagsPro->Fill(0.5,fCalculateSymmetryPlanes);
3035  //fSymmetryPlanesFlagsPro->GetXaxis()->SetBinLabel(2,"TBI"); fSymmetryPlanesFlagsPro->Fill(1.5,TBI);
3037 
3038  if(!fCalculateSymmetryPlanes){return;}
3039 
3040  // b) Book TProfile *fSymmetryPlanesPro[?][?]: TBI check the exact dimensions in the header file
3041  TString sTitle[1][2] = {"#LT#LTcos[4(#psi_{2}-#psi_{1})]#GT#GT","#LT#LTcos[4(#psi_{4}-#psi_{2})]#GT#GT"}; // TBI check h.w. harmonic 4, as well as h.w. indices on symmetry planes
3042  TString sLabels[4] = {"k = 0","k = 1","k = 2","k = 3"};
3043  for(Int_t gc=0;gc<1;gc++) // 'generic correlator': [[0]:(Psi2n,Psi1n),[1]:...] TBI upper boundary will change
3044  {
3045  for(Int_t n=0;n<2;n++) // 'harmonic n for generic correlator': [[0]:n=1,[1]:n=2,...] TBI upper boundary will change
3046  {
3047  fSymmetryPlanesPro[gc][n] = new TProfile(Form("%d,%d",gc,n),sTitle[gc][n].Data(),4,0.,4.); // TBI this is a landmine, introduce fnHighestOptimizerSPC eventually instead of '4'
3048  fSymmetryPlanesPro[gc][n]->Sumw2();
3049  fSymmetryPlanesPro[gc][n]->SetStats(kFALSE);
3050  fSymmetryPlanesPro[gc][n]->SetMarkerColor(kRed);
3051  fSymmetryPlanesPro[gc][n]->SetMarkerStyle(kFullSquare);
3052  fSymmetryPlanesPro[gc][n]->SetLineColor(kBlack);
3053  fSymmetryPlanesPro[gc][n]->SetLabelSize(0.0544);
3054  for(Int_t o=0;o<4;o++) // TBI h.w. '4'
3055  {
3056  fSymmetryPlanesPro[gc][n]->GetXaxis()->SetBinLabel(o+1,sLabels[o].Data());
3057  } // for(Int_t o=0;o<4;o++)
3059  }
3060  }
3061 
3062 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForSymmetryPlanes()
3063 
3064 //=======================================================================================================================
3065 
3067 {
3068  // Book all the stuff for correlations with eta gaps.
3069 
3070  // a) Book the profile holding all the flags for correlations with eta gaps;
3071  // b) Book TProfile *fEtaGapsPro[6][10]
3072 
3073  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForEtaGaps()";
3074 
3075  // a) Book the profile holding all the flags for correlations with eta gaps:
3076  fEtaGapsFlagsPro = new TProfile("fEtaGapsFlagsPro","Flags for correlations with eta gaps",1,0,1);
3077  fEtaGapsFlagsPro->SetTickLength(-0.01,"Y");
3078  fEtaGapsFlagsPro->SetMarkerStyle(25);
3079  fEtaGapsFlagsPro->SetLabelSize(0.03);
3080  fEtaGapsFlagsPro->SetLabelOffset(0.02,"Y");
3081  fEtaGapsFlagsPro->SetStats(kFALSE);
3082  fEtaGapsFlagsPro->SetFillColor(kGray);
3083  fEtaGapsFlagsPro->SetLineColor(kBlack);
3084  fEtaGapsFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateEtaGaps"); fEtaGapsFlagsPro->Fill(0.5,fCalculateEtaGaps);
3085  fEtaGapsFlagsPro->GetXaxis()->SetBinLabel(2,"fLowestHarmonicEtaGaps"); fEtaGapsFlagsPro->Fill(1.5,fLowestHarmonicEtaGaps);
3086  fEtaGapsFlagsPro->GetXaxis()->SetBinLabel(3,"fHighestHarmonicEtaGaps"); fEtaGapsFlagsPro->Fill(2.5,fHighestHarmonicEtaGaps);
3088 
3089  if(!fCalculateEtaGaps){return;}
3090 
3091  // b) Book TProfile *fEtaGapsPro[6][10]:
3092  Int_t markerColor[6] = {kBlue,kRed,kGreen+2,kBlack,kBlack,kBlack};
3093  TString sEtaGaps[11] = {"|#Delta#eta| > 1.0","|#Delta#eta| > 0.9","|#Delta#eta| > 0.8","|#Delta#eta| > 0.7","|#Delta#eta| > 0.6","|#Delta#eta| > 0.5","|#Delta#eta| > 0.4","|#Delta#eta| > 0.3","|#Delta#eta| > 0.2","|#Delta#eta| > 0.1","|#Delta#eta| > 0.0"};
3094  for(Int_t h=fLowestHarmonicEtaGaps-1;h<=fHighestHarmonicEtaGaps-1;h++) // harmonics
3095  {
3096  fEtaGapsPro[h] = new TProfile(Form("EtaGaps_v%d",h+1),Form("2-p correlation for harmonic v_{%d} with eta gaps",h+1),11,0.,11.);
3097  fEtaGapsPro[h]->Sumw2();
3098  fEtaGapsPro[h]->SetStats(kFALSE);
3099  fEtaGapsPro[h]->SetMarkerColor(markerColor[h]);
3100  fEtaGapsPro[h]->SetMarkerStyle(kFullSquare);
3101  fEtaGapsPro[h]->SetLineColor(markerColor[h]);
3102  fEtaGapsPro[h]->SetLabelSize(0.0544);
3103  fEtaGapsPro[h]->SetMinimum(0.0);
3104  for(Int_t eg=0;eg<11;eg++)
3105  {
3106  fEtaGapsPro[h]->GetXaxis()->SetBinLabel(eg+1,sEtaGaps[eg].Data());
3107  }
3108  fEtaGapsList->Add(fEtaGapsPro[h]);
3109  }
3110 
3111 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForEtaGaps()
3112 
3113 //=======================================================================================================================
3114 
3116 {
3117  // Book all the stuff for event-by-event cumulants.
3118 
3119  // a) Book the profile holding all the flags for e-b-e cumulants;
3120  // b) Book TProfile *fEbECumulantsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]).
3121 
3122  TString sMethodName = "void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForEbECumulants()";
3123 
3124  // a) Book the profile holding all the flags for e-b-e cumulants:
3125  fEbECumulantsFlagsPro = new TProfile("fEbECumulantsFlagsPro","Flags for e-b-e cumulants",1,0,1);
3126  fEbECumulantsFlagsPro->SetTickLength(-0.01,"Y");
3127  fEbECumulantsFlagsPro->SetMarkerStyle(25);
3128  fEbECumulantsFlagsPro->SetLabelSize(0.03);
3129  fEbECumulantsFlagsPro->SetLabelOffset(0.02,"Y");
3130  fEbECumulantsFlagsPro->SetStats(kFALSE);
3131  fEbECumulantsFlagsPro->SetFillColor(kGray);
3132  fEbECumulantsFlagsPro->SetLineColor(kBlack);
3133  fEbECumulantsFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateEbECumulants"); fEbECumulantsFlagsPro->Fill(0.5,fCalculateEbECumulants);
3135 
3136  if(!fCalculateEbECumulants){return;} // TBI is this safe enough?
3137 
3138  // b) Book TProfile *fEbECumulantsPro[2][8] ([0=cos,1=sin][1p,2p,...,8p]):
3139  TString sCosSin[2] = {"Cos","Sin"};
3140  Int_t markerColor[2] = {kBlue,kRed};
3141  Int_t markerStyle[2] = {kFullSquare,kFullSquare};
3142  Int_t nBins[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
3143  Int_t nBinsTitle[8] = {1,1,1,1,1,1,1,1}; // TBI hardwired 8, shall be fMaxCorrelator
3144  Int_t nToBeFilled[8] = {0,0,0,0,0,0,0,0}; // TBI hardwired 8, shall be fMaxCorrelator
3145  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
3146  {
3147  // Implementing \binom{n+k-1}{k}, which is the resulting number of sets obtained
3148  // after sampling n starting elements into k subsets, repetitions allowed.
3149  // In my case, n=2*fMaxHarmonic+1, k=c+1, hence:
3150  nBins[c] = (Int_t)(TMath::Factorial(2*fMaxHarmonic+1+c+1-1)
3151  / (TMath::Factorial(2*fMaxHarmonic+1-1)*TMath::Factorial(c+1)));
3152  nBinsTitle[c] = nBins[c];
3153  if(c>=fDontGoBeyond){nBins[c]=1;} // TBI a bit of spaghetti here...
3154  } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
3155  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
3156  {
3157  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
3158  {
3159  fEbECumulantsPro[cs][c] = new TProfile(Form("%dpEbECumulants%s",c+1,sCosSin[cs].Data()),"",nBins[c],0.,1.*nBins[c]);
3160  fEbECumulantsPro[cs][c]->Sumw2();
3161  fEbECumulantsPro[cs][c]->SetStats(kFALSE);
3162  fEbECumulantsPro[cs][c]->SetMarkerColor(markerColor[cs]);
3163  fEbECumulantsPro[cs][c]->SetMarkerStyle(markerStyle[cs]);
3164  fEbECumulantsList->Add(fEbECumulantsPro[cs][c]);
3165  } // for(Int_t c=0;c<8;c++) // [1p,2p,...,8p]
3166  } // for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
3167  // Set all bin labels: TBI this can be implemented better, most likely...
3168  Int_t binNo[8]; for(Int_t c=0;c<fMaxCorrelator;c++){binNo[c]=1;} // TBI hardwired 8, shall be fMaxCorrelator
3169  for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++)
3170  {
3171  if(fSkipZeroHarmonics && 0==n1){continue;}
3172  if(fCalculateAll)
3173  {
3174  fEbECumulantsPro[0][0]->GetXaxis()->SetBinLabel(binNo[0],Form("Cos(%d)",n1));
3175  fEbECumulantsPro[1][0]->GetXaxis()->SetBinLabel(binNo[0]++,Form("Sin(%d)",n1));
3176  nToBeFilled[0]++;
3177  }
3178  if(1==fDontGoBeyond){continue;}
3179  for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
3180  {
3181  if(fSkipZeroHarmonics && 0==n2){continue;}
3182  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2))
3183  || (fCalculateSameIsotropic && 0==n1+n2 && TMath::Abs(n1)==TMath::Abs(n2)))
3184  {
3185  fEbECumulantsPro[0][1]->GetXaxis()->SetBinLabel(binNo[1],Form("Cos(%d,%d)",n1,n2));
3186  fEbECumulantsPro[1][1]->GetXaxis()->SetBinLabel(binNo[1]++,Form("Sin(%d,%d)",n1,n2));
3187  nToBeFilled[1]++;
3188  }
3189  if(2==fDontGoBeyond){continue;}
3190  for(Int_t n3=n2;n3<=fMaxHarmonic;n3++)
3191  {
3192  if(fSkipZeroHarmonics && 0==n3){continue;}
3193  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3) || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3))
3194  || (fCalculateSameIsotropic && 0==n1+n2+n3 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)))
3195  {
3196  fEbECumulantsPro[0][2]->GetXaxis()->SetBinLabel(binNo[2],Form("Cos(%d,%d,%d)",n1,n2,n3));
3197  fEbECumulantsPro[1][2]->GetXaxis()->SetBinLabel(binNo[2]++,Form("Sin(%d,%d,%d)",n1,n2,n3));
3198  nToBeFilled[2]++;
3199  }
3200  if(3==fDontGoBeyond){continue;}
3201  for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)
3202  {
3203  if(fSkipZeroHarmonics && 0==n4){continue;}
3204  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))
3205  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)))
3206  {
3207  fEbECumulantsPro[0][3]->GetXaxis()->SetBinLabel(binNo[3],Form("Cos(%d,%d,%d,%d)",n1,n2,n3,n4));
3208  fEbECumulantsPro[1][3]->GetXaxis()->SetBinLabel(binNo[3]++,Form("Sin(%d,%d,%d,%d)",n1,n2,n3,n4));
3209  nToBeFilled[3]++;
3210  }
3211  if(4==fDontGoBeyond){continue;}
3212  for(Int_t n5=n4;n5<=fMaxHarmonic;n5++)
3213  {
3214  if(fSkipZeroHarmonics && 0==n5){continue;}
3215  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5)
3216  || (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))
3217  || (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)))
3218  {
3219  fEbECumulantsPro[0][4]->GetXaxis()->SetBinLabel(binNo[4],Form("Cos(%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5));
3220  fEbECumulantsPro[1][4]->GetXaxis()->SetBinLabel(binNo[4]++,Form("Sin(%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5));
3221  nToBeFilled[4]++;
3222  }
3223  if(5==fDontGoBeyond){continue;}
3224  for(Int_t n6=n5;n6<=fMaxHarmonic;n6++)
3225  {
3226  if(fSkipZeroHarmonics && 0==n6){continue;}
3227  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6)
3228  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
3229  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6))
3230  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
3231  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6)))
3232  {
3233  fEbECumulantsPro[0][5]->GetXaxis()->SetBinLabel(binNo[5],Form("Cos(%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6));
3234  fEbECumulantsPro[1][5]->GetXaxis()->SetBinLabel(binNo[5]++,Form("Sin(%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6));
3235  nToBeFilled[5]++;
3236  }
3237  if(6==fDontGoBeyond){continue;}
3238  for(Int_t n7=n6;n7<=fMaxHarmonic;n7++)
3239  {
3240  if(fSkipZeroHarmonics && 0==n7){continue;}
3241  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7)
3242  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
3243  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7))
3244  || (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)
3245  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7)))
3246  {
3247  fEbECumulantsPro[0][6]->GetXaxis()->SetBinLabel(binNo[6],Form("Cos(%d,%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6,n7));
3248  fEbECumulantsPro[1][6]->GetXaxis()->SetBinLabel(binNo[6]++,Form("Sin(%d,%d,%d,%d,%d,%d,%d)",n1,n2,n3,n4,n5,n6,n7));
3249  nToBeFilled[6]++;
3250  }
3251  if(7==fDontGoBeyond){continue;}
3252  for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
3253  {
3254  if(fSkipZeroHarmonics && 0==n8){continue;}
3255  if(fCalculateAll || (fCalculateIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8)
3256  || (fCalculateSame && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3) && TMath::Abs(n1)==TMath::Abs(n4)
3257  && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7) && TMath::Abs(n1)==TMath::Abs(n8))
3258  || (fCalculateSameIsotropic && 0==n1+n2+n3+n4+n5+n6+n7+n8 && TMath::Abs(n1)==TMath::Abs(n2) && TMath::Abs(n1)==TMath::Abs(n3)
3259  && TMath::Abs(n1)==TMath::Abs(n4) && TMath::Abs(n1)==TMath::Abs(n5) && TMath::Abs(n1)==TMath::Abs(n6) && TMath::Abs(n1)==TMath::Abs(n7)
3260  && TMath::Abs(n1)==TMath::Abs(n8)))
3261  {
3262  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));
3263  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));
3264  nToBeFilled[7]++;
3265  }
3266  } // for(Int_t n8=n7;n8<=fMaxHarmonic;n8++)
3267  } // for(Int_t n7=n6;n7<=fMaxHarmonic;n7++)
3268  } // for(Int_t n6=n5;n6<=fMaxHarmonic;n6++)
3269  } // for(Int_t n5=n4;n5<=fMaxHarmonic;n5++)
3270  } // for(Int_t n4=n3;n4<=fMaxHarmonic;n4++)
3271  } // for(Int_t n3=n2;n3<=fMaxHarmonic;n3++)
3272  } // for(Int_t n2=n1;n2<=fMaxHarmonic;n2++)
3273  } // for(Int_t n1=-fMaxHarmonic;n1<=fMaxHarmonic;n1++)
3274 
3275  for(Int_t cs=0;cs<2;cs++) // [0=cos,1=sin]
3276  {
3277  for(Int_t c=0;c<fMaxCorrelator;c++) // [1p,2p,...,8p]
3278  {
3279  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]));
3280  fEbECumulantsPro[cs][c]->GetXaxis()->SetRangeUser(0.,fEbECumulantsPro[cs][c]->GetBinLowEdge(nToBeFilled[c]+1));
3281  }
3282  }
3283 
3284 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForEbECumulants()
3285 
3286 //=======================================================================================================================
3287 
3289 {
3290  // Book all the stuff for nested loops.
3291 
3292  // TBI this method is just ugly, who implemented it like this...
3293 
3294  // a) Set default harmonic values;
3295  // b) Book the profile holding all the flags for nested loops;
3296  // c) Book the profile holding all results for nested loops (cosine);
3297  // d) Book the profile holding all results for nested loops (sinus);
3298  // e) Book the profile holding all results for differential nested loops.
3299 
3300  // a) Set default harmonic values:
3301  //delete gRandom; // TBI this is not safe here,
3302  //gRandom = new TRandom3(0);
3303  Int_t h1 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1); // TBI reimplement all these lines eventually
3304  Int_t h2 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
3305  Int_t h3 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
3306  Int_t h4 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
3307  Int_t h5 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
3308  Int_t h6 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
3309  Int_t h7 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
3310  Int_t h8 = (Int_t)pow(-1.,1.*gRandom->Integer(fMaxHarmonic+1))*gRandom->Integer(fMaxHarmonic+1);
3311 
3312  // REMARK: This values can be overriden in a steering macro via
3313  // mpc->GetNestedLoopsFlagsPro()->SetBinContent(<binNo>,<value>);
3314 
3315  // b) Book the profile holding all the flags for nested loops:
3316  fNestedLoopsFlagsPro = new TProfile("fNestedLoopsFlagsPro","Flags for nested loops",10,0,10);
3317  fNestedLoopsFlagsPro->SetTickLength(-0.01,"Y");
3318  fNestedLoopsFlagsPro->SetMarkerStyle(25);
3319  fNestedLoopsFlagsPro->SetLabelSize(0.03);
3320  fNestedLoopsFlagsPro->SetLabelOffset(0.02,"Y");
3321  fNestedLoopsFlagsPro->SetStats(kFALSE);
3322  fNestedLoopsFlagsPro->SetFillColor(kGray);
3323  fNestedLoopsFlagsPro->SetLineColor(kBlack);
3324  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(1,"fCrossCheckWithNestedLoops"); fNestedLoopsFlagsPro->Fill(0.5,fCrossCheckWithNestedLoops);
3325  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(2,"h_{1}"); fNestedLoopsFlagsPro->Fill(1.5,h1);
3326  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(3,"h_{2}"); fNestedLoopsFlagsPro->Fill(2.5,h2);
3327  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(4,"h_{3}"); fNestedLoopsFlagsPro->Fill(3.5,h3);
3328  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(5,"h_{4}"); fNestedLoopsFlagsPro->Fill(4.5,h4);
3329  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(6,"h_{5}"); fNestedLoopsFlagsPro->Fill(5.5,h5);
3330  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(7,"h_{6}"); fNestedLoopsFlagsPro->Fill(6.5,h6);
3331  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(8,"h_{7}"); fNestedLoopsFlagsPro->Fill(7.5,h7);
3332  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(9,"h_{8}"); fNestedLoopsFlagsPro->Fill(8.5,h8);
3333  fNestedLoopsFlagsPro->GetXaxis()->SetBinLabel(10,"fCrossCheckDiffWithNestedLoops"); fNestedLoopsFlagsPro->Fill(9.5,fCrossCheckDiffWithNestedLoops);
3335 
3336  // c) Book the profile holding all results for nested loops (cosine):
3337  fNestedLoopsResultsCosPro = new TProfile("fNestedLoopsResultsCosPro","Nested loops results (cosine)",16,0,16);
3338  fNestedLoopsResultsCosPro->SetTickLength(-0.01,"Y");
3339  fNestedLoopsResultsCosPro->SetMarkerStyle(25);
3340  fNestedLoopsResultsCosPro->SetLabelSize(0.02);
3341  fNestedLoopsResultsCosPro->SetLabelOffset(0.02,"Y");
3342  fNestedLoopsResultsCosPro->SetStats(kFALSE);
3343  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(1,Form("N: 1p(%d)",h1));
3344  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(2,Form("Q: 1p(%d)",h1));
3345  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(3,Form("N: 2p(%d,%d)",h1,h2));
3346  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(4,Form("Q: 2p(%d,%d)",h1,h2));
3347  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(5,Form("N: 3p(%d,%d,%d)",h1,h2,h3));
3348  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(6,Form("Q: 3p(%d,%d,%d)",h1,h2,h3));
3349  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(7,Form("N: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
3350  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(8,Form("Q: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
3351  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(9,Form("N: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
3352  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(10,Form("Q: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
3353  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(11,Form("N: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
3354  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(12,Form("Q: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
3355  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(13,Form("N: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
3356  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(14,Form("Q: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
3357  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(15,Form("N: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
3358  fNestedLoopsResultsCosPro->GetXaxis()->SetBinLabel(16,Form("Q: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
3360 
3361  // d) Book the profile holding all results for nested loops (sinus):
3362  fNestedLoopsResultsSinPro = new TProfile("fNestedLoopsResultsSinPro","Nested loops results (sinus)",16,0,16);
3363  fNestedLoopsResultsSinPro->SetTickLength(-0.01,"Y");
3364  fNestedLoopsResultsSinPro->SetMarkerStyle(25);
3365  fNestedLoopsResultsSinPro->SetLabelSize(0.02);
3366  fNestedLoopsResultsSinPro->SetLabelOffset(0.02,"Y");
3367  fNestedLoopsResultsSinPro->SetStats(kFALSE);
3368  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(1,Form("N: 1p(%d)",h1));
3369  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(2,Form("Q: 1p(%d)",h1));
3370  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(3,Form("N: 2p(%d,%d)",h1,h2));
3371  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(4,Form("Q: 2p(%d,%d)",h1,h2));
3372  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(5,Form("N: 3p(%d,%d,%d)",h1,h2,h3));
3373  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(6,Form("Q: 3p(%d,%d,%d)",h1,h2,h3));
3374  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(7,Form("N: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
3375  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(8,Form("Q: 4p(%d,%d,%d,%d)",h1,h2,h3,h4));
3376  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(9,Form("N: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
3377  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(10,Form("Q: 5p(%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5));
3378  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(11,Form("N: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
3379  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(12,Form("Q: 6p(%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6));
3380  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(13,Form("N: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
3381  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(14,Form("Q: 7p(%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7));
3382  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(15,Form("N: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
3383  fNestedLoopsResultsSinPro->GetXaxis()->SetBinLabel(16,Form("Q: 8p(%d,%d,%d,%d,%d,%d,%d,%d)",h1,h2,h3,h4,h5,h6,h7,h8));
3385 
3386  // e) Book the profile holding all results for differential nested loops:
3387  fNestedLoopsDiffResultsPro = new TProfile("fNestedLoopsDiffResultsPro","Differential nested loops results",1,0.,1.);
3388  fNestedLoopsDiffResultsPro->SetStats(kFALSE);
3390 
3391 } // void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForNestedLoops()
3392 
3393 //=======================================================================================================================
3394 
3396 {
3397  // Book all the stuff for 'standard candles'.
3398 
3399  // a) Book the profile holding all the flags for 'standard candles';
3400  // b) Book the histogram holding all results for 'standard candles';
3401  // c) Book TProfile2D *fProductsSCPro.
3402 
3403  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForStandardCandles()";
3404 
3405  // a) Book the profile holding all the flags for 'standard candles':
3406  fStandardCandlesFlagsPro = new TProfile("fStandardCandlesFlagsPro","Flags for standard candles",2,0,2);
3407  fStandardCandlesFlagsPro->SetTickLength(-0.01,"Y");
3408  fStandardCandlesFlagsPro->SetMarkerStyle(25);
3409  fStandardCandlesFlagsPro->SetLabelSize(0.03);
3410  fStandardCandlesFlagsPro->SetLabelOffset(0.02,"Y");
3411  fStandardCandlesFlagsPro->SetStats(kFALSE);
3412  fStandardCandlesFlagsPro->SetFillColor(kGray);
3413  fStandardCandlesFlagsPro->SetLineColor(kBlack);
3414  fStandardCandlesFlagsPro->GetXaxis()->SetBinLabel(1,"fCalculateStandardCandles"); fStandardCandlesFlagsPro->Fill(0.5,fCalculateStandardCandles);
3415  fStandardCandlesFlagsPro->GetXaxis()->SetBinLabel(2,"fPropagateErrorSC"); fStandardCandlesFlagsPro->Fill(1.5,fPropagateErrorSC);
3417 
3418  if(!fCalculateStandardCandles){return;} // TBI is this safe like this?
3419 
3420  // b) Book the histogram holding all results for 'standard candles':
3421  Int_t nBins = fMaxHarmonic*(fMaxHarmonic-1)/2;
3422  fStandardCandlesHist = new TH1D("fStandardCandlesHist","Standard candles (SC)",nBins,0.,1.*nBins);
3423  fStandardCandlesHist->SetStats(kFALSE);
3424  fStandardCandlesHist->SetMarkerStyle(kFullSquare);
3425  fStandardCandlesHist->SetMarkerColor(kBlue);
3426  fStandardCandlesHist->SetLineColor(kBlue);
3427  Int_t binNo = 1;
3428  for(Int_t n1=-fMaxHarmonic;n1<=-2;n1++)
3429  {
3430  for(Int_t n2=n1+1;n2<=-1;n2++)
3431  {
3432  fStandardCandlesHist->GetXaxis()->SetBinLabel(binNo++,Form("SC(%d,%d,%d,%d)",n1,n2,-1*n2,-1*n1));
3433  }
3434  }
3435  if(binNo-1 != nBins){Fatal(sMethodName.Data(),"Well, binNo-1 != nBins ... :'( ");}
3436  fStandardCandlesList->Add(fStandardCandlesHist);
3437 
3438  if(!fPropagateErrorSC){return;} // TBI is this safe like this?
3439 
3440  // c) Book TProfile2D *fProductsSCPro:
3441  Int_t nBins2D = fMaxHarmonic + fMaxHarmonic*(fMaxHarmonic-1)/2; // #2-p + #4-p distinct correlators in SC context
3442  if(nBins2D<=0){Fatal(sMethodName.Data(),"nBins2D<=0");} // well, just in case...
3443  fProductsSCPro = new TProfile2D("fProductsSCPro","Products of correlations",nBins2D,0.,nBins2D,nBins2D,0.,nBins2D);
3444  fProductsSCPro->SetStats(kFALSE);
3445  fProductsSCPro->Sumw2();
3446  for(Int_t b=1;b<=fMaxHarmonic;b++) // 2-p correlators
3447  {
3448  fProductsSCPro->GetXaxis()->SetBinLabel(b,Form("Cos(%d,%d)",-(fMaxHarmonic+1-b),(fMaxHarmonic+1-b)));
3449  fProductsSCPro->GetYaxis()->SetBinLabel(b,Form("Cos(%d,%d)",-(fMaxHarmonic+1-b),(fMaxHarmonic+1-b)));
3450  } // for(Int_t b=1;b<=fMaxHarmonic;b++) // 2-p correlators
3451  for(Int_t b=fMaxHarmonic+1;b<=nBins2D;b++) // 4-p correlators
3452  {
3453  TString sBinLabel = fStandardCandlesHist->GetXaxis()->GetBinLabel(b-fMaxHarmonic);
3454  if(sBinLabel.EqualTo("")){Fatal(sMethodName.Data(),"sBinLabel.EqualTo...");}
3455  fProductsSCPro->GetXaxis()->SetBinLabel(b,sBinLabel.ReplaceAll("SC","Cos").Data());
3456  fProductsSCPro->GetYaxis()->SetBinLabel(b,sBinLabel.ReplaceAll("SC","Cos").Data());
3457  } // for(Int_t b=fMaxHarmonic+1;b<=nBins2D;b++) // 4-p correlators
3458  fStandardCandlesList->Add(fProductsSCPro);
3459 
3460 } // end of void AliFlowAnalysisWithMultiparticleCorrelations::BookEverythingForStandardCandles()
3461 
3462 //=======================================================================================================================
3463 
3465 {
3466  // Get pointers for everything and everywhere from the base list "fHistList".
3467 
3468  // a) Get pointer for base list fHistList;
3469  // b) Get pointer for profile holding internal flags and, well, set again all flags;
3470  // c) Get pointers for control histograms;
3471  // d) Get pointers for Q-vector;
3472  // e) Get pointers for correlations;
3473  // f) Get pointers for 'standard candles';
3474  // g) Get pointers for Q-cumulants;
3475  // h) Get pointers for differential correlations;
3476  // i) Get pointers for symmetry planes.
3477 
3478  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetOutputHistograms(TList *histList)";
3479 
3480  // a) Get pointer for base list fHistList and profile holding internal flags;
3481  fHistList = histList;
3482  if(!fHistList){Fatal(sMethodName.Data(),"fHistList is malicious today...");}
3483 
3484  // b) Get pointer for profile holding internal flags and, well, set again all flags:
3485  fInternalFlagsPro = dynamic_cast<TProfile*>(fHistList->FindObject("fInternalFlagsPro"));
3486  if(!fInternalFlagsPro){Fatal(sMethodName.Data(),"fInternalFlagsPro");}
3487  fUseInternalFlags = fInternalFlagsPro->GetBinContent(1);
3488  fMinNoRPs = (Int_t)fInternalFlagsPro->GetBinContent(2);
3489  fMaxNoRPs = (Int_t)fInternalFlagsPro->GetBinContent(3);
3490  fExactNoRPs = (Int_t)fInternalFlagsPro->GetBinContent(4);
3491  fPropagateError = (Bool_t)fInternalFlagsPro->GetBinContent(5);
3492 
3493  // c) Get pointers for control histograms:
3495 
3496  // d) Get pointers for Q-vector:
3497  this->GetPointersForQvector();
3498 
3499  // e) Get pointers for correlations:
3500  this->GetPointersForCorrelations();
3501 
3502  // f) Get pointers for 'standard candles':
3504 
3505  // g) Get pointers for Q-cumulants:
3506  this->GetPointersForQcumulants();
3507 
3508  // h) Get pointers for differential correlations:
3510 
3511  // i) Get pointers for symmetry planes:
3513 
3514 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetOutputHistograms(TList *histList)
3515 
3516 //=======================================================================================================================
3517 
3519 {
3520  // Get pointers for Q-vector objects.
3521 
3522  // a) Get pointer for fQvectorList;
3523  // b) Get pointer for fQvectorFlagsPro;
3524  // c) Set again all flags.
3525 
3526  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForQvector()";
3527 
3528  // a) Get pointer for fQvectorList:
3529  fQvectorList = dynamic_cast<TList*>(fHistList->FindObject("Q-vectors"));
3530  if(!fQvectorList){Fatal(sMethodName.Data(),"fQvectorList");}
3531 
3532  // b) Get pointer for fQvectorFlagsPro:
3533  fQvectorFlagsPro = dynamic_cast<TProfile*>(fQvectorList->FindObject("fQvectorFlagsPro"));
3534  if(!fQvectorFlagsPro){Fatal(sMethodName.Data(),"fQvectorFlagsPro");}
3535 
3536  // c) Set again all flags:
3537  fCalculateQvector = (Bool_t)fQvectorFlagsPro->GetBinContent(1);
3538  fCalculateDiffQvectors = (Bool_t)fQvectorFlagsPro->GetBinContent(2);
3539 
3540 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForQvector()
3541 
3542 //=======================================================================================================================
3543 
3545 {
3546  // Get pointers for 'standard candles'.
3547 
3548  // a) Get pointer for fStandardCandlesList;
3549  // b) Get pointer for fStandardCandlesFlagsPro;
3550  // c) Set again all flags;
3551  // d) Get pointer TH1D *fStandardCandlesHist;
3552  // e) Get pointer TProfile2D *fProductsSCPro.
3553 
3554  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForStandardCandles()";
3555 
3556  // a) Get pointer for fStandardCandlesList:
3557  fStandardCandlesList = dynamic_cast<TList*>(fHistList->FindObject("Standard Candles"));
3558  if(!fStandardCandlesList){Fatal(sMethodName.Data(),"fStandardCandlesList");}
3559 
3560  // b) Get pointer for fStandardCandlesFlagsPro:
3561  fStandardCandlesFlagsPro = dynamic_cast<TProfile*>(fStandardCandlesList->FindObject("fStandardCandlesFlagsPro"));
3562  if(!fStandardCandlesFlagsPro){Fatal(sMethodName.Data(),"fStandardCandlesFlagsPro");}
3563 
3564  // c) Set again all flags:
3565  fCalculateStandardCandles = (Bool_t)fStandardCandlesFlagsPro->GetBinContent(1);
3566  fPropagateErrorSC = (Bool_t)fStandardCandlesFlagsPro->GetBinContent(2);
3567 
3568  if(!fCalculateStandardCandles){return;} // TBI is this safe enough
3569 
3570  // d) Get pointer TH1D *fStandardCandlesHist:
3571  fStandardCandlesHist = dynamic_cast<TH1D*>(fStandardCandlesList->FindObject("fStandardCandlesHist"));
3572 
3573  if(!fPropagateErrorSC){return;} // TBI is this safe enough
3574 
3575  // e) Get pointer TProfile2D *fProductsSCPro:
3576  fProductsSCPro = dynamic_cast<TProfile2D*>(fStandardCandlesList->FindObject("fProductsSCPro"));
3577 
3578 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForStandardCandles()
3579 
3580 //=======================================================================================================================
3581 
3583 {
3584  // Get pointers for control histograms.
3585 
3586  // a) Get pointer for fControlHistogramsList; TBI
3587  // b) Get pointer for fControlHistogramsFlagsPro; TBI
3588  // c) Set again all flags; TBI
3589  // d) Get pointers to TH1D *fKinematicsHist[2][3]; TBI
3590  // e) Get pointers to TH1D *fMultDistributionsHist[3]; TBI
3591  // f) Get pointers to TH2D *fMultCorrelationsHist[3]. TBI
3592 
3593  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForControlHistograms()";
3594 
3595  // a) Get pointer for fControlHistogramsList: TBI
3596  fControlHistogramsList = dynamic_cast<TList*>(fHistList->FindObject("Control Histograms"));
3597  if(!fControlHistogramsList){Fatal(sMethodName.Data(),"fControlHistogramsList");}
3598 
3599  // b) Get pointer for fControlHistogramsFlagsPro: TBI
3600  fControlHistogramsFlagsPro = dynamic_cast<TProfile*>(fControlHistogramsList->FindObject("fControlHistogramsFlagsPro"));
3601  if(!fControlHistogramsFlagsPro){Fatal(sMethodName.Data(),"fControlHistogramsFlagsPro");}
3602 
3603  // c) Set again all flags:
3604  fFillControlHistograms = fControlHistogramsFlagsPro->GetBinContent(1);
3605  fFillKinematicsHist = fControlHistogramsFlagsPro->GetBinContent(2);
3606  fFillMultDistributionsHist = fControlHistogramsFlagsPro->GetBinContent(3);
3607  fFillMultCorrelationsHist = fControlHistogramsFlagsPro->GetBinContent(4);
3608  fDontFill[0] = fControlHistogramsFlagsPro->GetBinContent(5);
3609  fDontFill[1] = fControlHistogramsFlagsPro->GetBinContent(6);
3610  fDontFill[2] = fControlHistogramsFlagsPro->GetBinContent(7);
3611 
3612  if(!fFillControlHistograms){return;} // TBI is this safe enough
3613 
3614  // d) Get pointers to fKinematicsHist[2][3]: TBI
3615  TString name[2][3] = {{"RP,phi","RP,pt","RP,eta"},{"POI,phi","POI,pt","POI,eta"}}; // [RP,POI][phi,pt,eta]
3616  for(Int_t rp=0;rp<2;rp++) // [RP,POI]
3617  {
3618  if(fDontFill[rp]){continue;}
3619  for(Int_t ppe=0;ppe<3;ppe++) // [phi,pt,eta]
3620  {
3621  fKinematicsHist[rp][ppe] = dynamic_cast<TH1D*>(fControlHistogramsList->FindObject(name[rp][ppe].Data()));
3622  if(!fKinematicsHist[rp][ppe] && fFillKinematicsHist){Fatal(sMethodName.Data(),"%s",name[rp][ppe].Data());} // TBI
3623  }
3624  }
3625 
3626  // e) Get pointers to TH1D *fMultDistributionsHist[3]:
3627  TString nameMult[3] = {"Multiplicity (RP)","Multiplicity (POI)","Multiplicity (REF)"}; // [RP,POI,reference multiplicity]
3628  for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
3629  {
3630  if(fDontFill[rprm]){continue;}
3631  fMultDistributionsHist[rprm] = dynamic_cast<TH1D*>(fControlHistogramsList->FindObject(nameMult[rprm].Data()));
3632  if(!fMultDistributionsHist[rprm] && fFillMultDistributionsHist){Fatal(sMethodName.Data(),"%s",nameMult[rprm].Data());} // TBI
3633  } // for(Int_t rprm=0;rprm<3;rprm++) // [RP,POI,reference multiplicity]
3634 
3635  // f) Get pointers to TH2I *fMultCorrelationsHist[3]: TBI automatize the things here (at some point...)...
3636  if(!fDontFill[0] && !fDontFill[1])
3637  {
3638  fMultCorrelationsHist[0] = dynamic_cast<TH2I*>(fControlHistogramsList->FindObject("Multiplicity (RP vs. POI)"));
3639  if(!fMultCorrelationsHist[0] && fFillMultCorrelationsHist){Fatal(sMethodName.Data(),"Multiplicity (RP vs. POI)");} // TBI
3640  }
3641  if(!fDontFill[0] && !fDontFill[2])
3642  {
3643  fMultCorrelationsHist[1] = dynamic_cast<TH2I*>(fControlHistogramsList->FindObject("Multiplicity (RP vs. REF)"));
3644  if(!fMultCorrelationsHist[1] && fFillMultCorrelationsHist){Fatal(sMethodName.Data(),"Multiplicity (RP vs. REF)");} // TBI
3645  }
3646  if(!fDontFill[1] && !fDontFill[2])
3647  {
3648  fMultCorrelationsHist[2] = dynamic_cast<TH2I*>(fControlHistogramsList->FindObject("Multiplicity (POI vs. REF)"));
3649  if(!fMultCorrelationsHist[2] && fFillMultCorrelationsHist){Fatal(sMethodName.Data(),"Multiplicity (POI vs. REF)");} // TBI
3650  }
3651 
3652 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForControlHistograms()
3653 
3654 //=======================================================================================================================
3655 
3657 {
3658  // Get pointers for correlations.
3659 
3660  // a) Get pointer for fCorrelationsList; TBI
3661  // b) Get pointer for fCorrelationsFlagsPro; TBI
3662  // c) Set again all flags; TBI
3663  // d) Get pointers to TProfile *fCorrelationsPro[2][8].
3664 
3665  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForCorrelations()";
3666 
3667  // a) Get pointer for fCorrelationsList: TBI
3668  fCorrelationsList = dynamic_cast<TList*>(fHistList->FindObject("Correlations"));
3669  if(!fCorrelationsList){Fatal(sMethodName.Data(),"fCorrelationsList");}
3670 
3671  // b) Get pointer for fCorrelationsFlagsPro: TBI
3672  fCorrelationsFlagsPro = dynamic_cast<TProfile*>(fCorrelationsList->FindObject("fCorrelationsFlagsPro"));
3673 
3674  if(!fCorrelationsFlagsPro){Fatal(sMethodName.Data(),"fCorrelationsFlagsPro");}
3675 
3676  // c) Set again all flags:
3677  fCalculateCorrelations = fCorrelationsFlagsPro->GetBinContent(1);
3678  fMaxHarmonic = (Int_t)fCorrelationsFlagsPro->GetBinContent(2);
3679  fMaxCorrelator = (Int_t)fCorrelationsFlagsPro->GetBinContent(3);
3680  fCalculateIsotropic = (Bool_t)fCorrelationsFlagsPro->GetBinContent(4);
3681  fCalculateSame = (Bool_t)fCorrelationsFlagsPro->GetBinContent(5);
3682  fSkipZeroHarmonics = (Bool_t)fCorrelationsFlagsPro->GetBinContent(6);
3683  fCalculateSameIsotropic = (Bool_t)fCorrelationsFlagsPro->GetBinContent(7);
3684  fCalculateAll = (Bool_t)fCorrelationsFlagsPro->GetBinContent(8);
3685  fDontGoBeyond = (Int_t)fCorrelationsFlagsPro->GetBinContent(9);
3686  fCalculateOnlyForHarmonicQC = (Bool_t)fCorrelationsFlagsPro->GetBinContent(10);
3687  fCalculateOnlyForSC = (Bool_t)fCorrelationsFlagsPro->GetBinContent(11);
3688  fCalculateOnlyCos = (Bool_t)fCorrelationsFlagsPro->GetBinContent(12);
3689  fCalculateOnlySin = (Bool_t)fCorrelationsFlagsPro->GetBinContent(13);
3690 
3691  if(!fCalculateCorrelations){return;} // TBI is this safe enough, that is the question...
3692 
3693  // d) Get pointers to TProfile *fCorrelationsPro[2][8]:
3694  TString sCosSin[2] = {"Cos","Sin"};
3695  for(Int_t cs=0;cs<2;cs++)
3696  {
3697  if(fCalculateOnlyCos && 1==cs){continue;}
3698  else if(fCalculateOnlySin && 0==cs){continue;}
3699  for(Int_t c=0;c<fMaxCorrelator;c++)
3700  {
3701  if(c==fDontGoBeyond){continue;}
3702  if(fCalculateOnlyForHarmonicQC && c%2==0){continue;}
3703  fCorrelationsPro[cs][c] = dynamic_cast<TProfile*>(fCorrelationsList->FindObject(Form("%dpCorrelations%s",c+1,sCosSin[cs].Data())));
3704  if(!fCorrelationsPro[cs][c]){Fatal(sMethodName.Data(),"%dpCorrelations%s",c+1,sCosSin[cs].Data());}
3705  }
3706  }
3707 
3708 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForCorrelations()
3709 
3710 //=======================================================================================================================
3711 
3713 {
3714  // Get pointers for Q-cumulants.
3715 
3716  // a) Get pointer for fQcumulantsList;
3717  // b) Get pointer for fQcumulantsFlagsPro;
3718  // c) Set again all flags;
3719  // d) Get pointer for fQcumulantsHist;
3720  // e) Get pointer for fReferenceFlowHist;
3721  // f) Get pointer for fProductsQCPro.
3722 
3723  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForQcumulants()";
3724 
3725  // a) Get pointer for fQcumulantsList:
3726  fQcumulantsList = dynamic_cast<TList*>(fHistList->FindObject("Q-cumulants"));
3727  if(!fQcumulantsList){Fatal(sMethodName.Data(),"fQcumulantsList");}
3728 
3729  // b) Get pointer for fQcumulantsFlagsPro:
3730  fQcumulantsFlagsPro = dynamic_cast<TProfile*>(fQcumulantsList->FindObject("fQcumulantsFlagsPro"));
3731  if(!fQcumulantsFlagsPro){Fatal(sMethodName.Data(),"fQcumulantsFlagsPro");}
3732 
3733  // c) Set again all flags:
3734  fCalculateQcumulants = (Bool_t) fQcumulantsFlagsPro->GetBinContent(1);
3735  fHarmonicQC = (Int_t) fQcumulantsFlagsPro->GetBinContent(2);
3736  fPropagateErrorQC = (Bool_t) fQcumulantsFlagsPro->GetBinContent(3);
3737 
3738  if(!fCalculateQcumulants){return;} // TBI is this safe enough
3739 
3740  // d) Get pointer for fQcumulantsHist:
3741  fQcumulantsHist = dynamic_cast<TH1D*>(fQcumulantsList->FindObject("Q-cumulants"));
3742  if(!fQcumulantsHist){Fatal(sMethodName.Data(),"fQcumulantsHist");}
3743 
3744  // e) Get pointer for fReferenceFlowHist:
3745  fReferenceFlowHist = dynamic_cast<TH1D*>(fQcumulantsList->FindObject("Reference Flow"));
3746  if(!fReferenceFlowHist){Fatal(sMethodName.Data(),"fReferenceFlowHist");}
3747 
3748  if(!fPropagateErrorQC){return;} // TBI is this safe enough
3749 
3750  // f) Get pointer for fProductsQCPro:
3751  fProductsQCPro = dynamic_cast<TProfile2D*>(fQcumulantsList->FindObject("fProductsQCPro"));
3752  if(!fProductsQCPro){Fatal(sMethodName.Data(),"fProductsQCPro");}
3753 
3754 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForQcumulants()
3755 
3756 //=======================================================================================================================
3757 
3759 {
3760  // Get pointers for differential correlations.
3761 
3762  // a) Get pointer for fDiffCorrelationsList;
3763  // b) Get pointer for fDiffCorrelationsFlagsPro;
3764  // c) Set again all flags;
3765  // d) Get pointers to TProfile *fDiffCorrelationsPro[2][4].
3766 
3767  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForDiffCorrelations()";
3768 
3769  // a) Get pointer for fDiffCorrelationsList:
3770  fDiffCorrelationsList = dynamic_cast<TList*>(fHistList->FindObject("Differential Correlations"));
3771  if(!fDiffCorrelationsList){Fatal(sMethodName.Data(),"fDiffCorrelationsList");}
3772 
3773  // b) Get pointer for fDiffCorrelationsFlagsPro:
3774  fDiffCorrelationsFlagsPro = dynamic_cast<TProfile*>(fDiffCorrelationsList->FindObject("fDiffCorrelationsFlagsPro"));
3775  if(!fDiffCorrelationsFlagsPro){Fatal(sMethodName.Data(),"fDiffCorrelationsFlagsPro");}
3776 
3777  // c) Set again all flags:
3778  fCalculateDiffCorrelations = fDiffCorrelationsFlagsPro->GetBinContent(1);
3779 
3780  if(!fCalculateDiffCorrelations){return;}
3781 
3782  // TBI get all pointers below for diff. stuff eventually, not needed for the time being.
3783 
3784  // d) Get pointers to TProfile *fDiffCorrelationsPro[2][4]: // TBI
3785  /*
3786  TString sCosSin[2] = {"Cos","Sin"};
3787  for(Int_t cs=0;cs<2;cs++)
3788  {
3789  if(fCalculateOnlyCos && 1==cs){continue;}
3790  else if(fCalculateOnlySin && 0==cs){continue;}
3791  for(Int_t c=0;c<fMaxCorrelator;c++)
3792  {
3793  fCorrelationsPro[cs][c] = dynamic_cast<TProfile*>(fCorrelationsList->FindObject(Form("%dpCorrelations%s",c+1,sCosSin[cs].Data())));
3794  if(!fCorrelationsPro[cs][c]){Fatal(sMethodName.Data(),"%dpCorrelations%s",c+1,sCosSin[cs].Data());}
3795  }
3796  }
3797  */
3798 
3799 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForDiffCorrelations()
3800 
3801 //=======================================================================================================================
3802 
3804 {
3805  // Get pointers for symmetry planes.
3806 
3807  // a) Get pointer for fSymmetryPlanesList;
3808  // b) Get pointer for fSymmetryPlanesFlagsPro;
3809  // c) Set again all flags;
3810  // d) Get pointers to TProfile *fSymmetryPlanesPro[?][?].
3811 
3812  TString sMethodName = "AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForSymmetryPlanes()";
3813 
3814  // a) Get pointer for fSymmetryPlanesList:
3815  fSymmetryPlanesList = dynamic_cast<TList*>(fHistList->FindObject("Symmetry_Plane_Correlations"));
3816  if(!fSymmetryPlanesList){Fatal(sMethodName.Data(),"fSymmetryPlanesList");}
3817 
3818  // b) Get pointer for fSymmetryPlanesFlagsPro:
3819  fSymmetryPlanesFlagsPro = dynamic_cast<TProfile*>(fSymmetryPlanesList->FindObject("fSymmetryPlanesFlagsPro"));
3820  if(!fSymmetryPlanesFlagsPro){Fatal(sMethodName.Data(),"fSymmetryPlanesFlagsPro");}
3821 
3822  // c) Set again all flags:
3823  fCalculateSymmetryPlanes = fSymmetryPlanesFlagsPro->GetBinContent(1);
3824 
3825  if(!fCalculateSymmetryPlanes){return;}
3826 
3827  // d) Get pointers to TProfile *fSymmetryPlanesPro[?][?]: // TBI
3828  for(Int_t gc=0;gc<1;gc++) // 'generic correlator': [[0]:(Psi2n,Psi1n),[1]:...] TBI upper boundary will change
3829  {
3830  for(Int_t n=0;n<2;n++) // 'harmonic n for generic correlator': [[0]:n=1,[1]:n=2,...] TBI upper boundary will change
3831  {
3832  fSymmetryPlanesPro[gc][n] = dynamic_cast<TProfile*>(fSymmetryPlanesList->FindObject(Form("%d,%d",gc,n)));
3833  if(!fSymmetryPlanesPro[gc][n]){Fatal(sMethodName.Data(),"fSymmetryPlanesPro[%d][%d]",gc,n);}
3834  }
3835  }
3836 
3837 } // void AliFlowAnalysisWithMultiparticleCorrelations::GetPointersForSymmetryPlanes()
3838 
3839 //=======================================================================================================================
3840 
3842 {
3843  // Initialize all arrays for Q-vector.
3844 
3845  for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
3846  {
3847  for(Int_t wp=0;wp<fMaxCorrelator+1;wp++) // weight power
3848  {
3849  fQvector[h][wp] = TComplex(0.,0.);
3850  for(Int_t b=0;b<100;b++) // TBI hardwired 100
3851  {
3852  fpvector[b][h][wp] = TComplex(0.,0.);
3853  fqvector[b][h][wp] = TComplex(0.,0.);
3854  }
3855  }
3856  }
3857 
3858 } // void AliFlowAnalysisWithMultiparticleCorrelations::InitializeArraysForQvector()
3859 
3860 //=======================================================================================================================
3861 
3863 {
3864  // Reset all Q-vector components to zero before starting a new event.
3865 
3866  for(Int_t h=0;h<fMaxHarmonic*fMaxCorrelator+1;h++)
3867  {
3868  for(Int_t wp=0;wp<fMaxCorrelator+1;wp++) // weight powe
3869  {
3870  fQvector[h][wp] = TComplex(0.,0.);
3871  if(!fCalculateDiffQvectors){continue;}
3872  for(Int_t b=0;b<100;b++) // TBI hardwired 100
3873  {
3874  fpvector[b][h][wp] = TComplex(0.,0.);
3875  fqvector[b][h][wp] = TComplex(0.,0.);
3876  }
3877  }
3878  }
3879 
3880 } // void AliFlowAnalysisWithMultiparticleCorrelations::ResetQvector()
3881 
3882 //=======================================================================================================================
3883 
3885 {
3886  // Using the fact that Q{-n,p} = Q{n,p}^*.
3887 
3888  if(n>=0){return fQvector[n][wp];}
3889  return TComplex::Conjugate(fQvector[-n][wp]);
3890 
3891 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::Q(Int_t n, Int_t wp)
3892 
3893 //=======================================================================================================================
3894 
3896 {
3897  // Using the fact that p{-n,p} = p{n,p}^*.
3898 
3899  if(n>=0){return fpvector[fDiffBinNo][n][wp];}
3900  return TComplex::Conjugate(fpvector[fDiffBinNo][-n][wp]);
3901 
3902 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::p(Int_t n, Int_t p)
3903 
3904 //=======================================================================================================================
3905 
3907 {
3908  // Using the fact that q{-n,p} = q{n,p}^*.
3909 
3910  // When weights are used for RPs and not for POIs, and vice versa, some gymnastics is required here:
3911  // TBI rethink and reimplement the lines below:
3912  Int_t nUseWeightsForRP = (Int_t)(fUseWeights[0][0] || fUseWeights[0][1] || fUseWeights[0][2]);
3913  Int_t nUseWeightsForPOI = (Int_t)(fUseWeights[1][0] || fUseWeights[1][1] || fUseWeights[1][2]);
3914  if(nUseWeightsForPOI == 1 && nUseWeightsForRP == 0){wp=1;}
3915  else if(nUseWeightsForPOI == 0 && nUseWeightsForRP == 1){wp-=1;}
3916 
3917  if(n>=0){return fqvector[fDiffBinNo][n][wp];}
3918  return TComplex::Conjugate(fqvector[fDiffBinNo][-n][wp]);
3919 
3920 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::q(Int_t n, Int_t wp)
3921 
3922 //=======================================================================================================================
3923 
3925 {
3926  // Generic expression <exp[i(n1*phi1)]>. TBI comment
3927 
3928  TComplex one = Q(n1,1);
3929 
3930  return one;
3931 
3932 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::One(Int_t n1)
3933 
3934 //=======================================================================================================================
3935 
3937 {
3938  // Generic two-particle correlation <exp[i(n1*phi1+n2*phi2)]>.
3939 
3940  TComplex two = Q(n1,1)*Q(n2,1)-Q(n1+n2,2);
3941 
3942  return two;
3943 
3944 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::Two(Int_t n1, Int_t n2)
3945 
3946 //=======================================================================================================================
3947 
3949 {
3950  // Generic three-particle correlation <exp[i(n1*phi1+n2*phi2+n3*phi3)]>.
3951 
3952  TComplex three = Q(n1,1)*Q(n2,1)*Q(n3,1)-Q(n1+n2,2)*Q(n3,1)-Q(n2,1)*Q(n1+n3,2)
3953  - Q(n1,1)*Q(n2+n3,2)+2.*Q(n1+n2+n3,3);
3954 
3955  return three;
3956 
3957 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::Three(Int_t n1, Int_t n2, Int_t n3)
3958 
3959 //=======================================================================================================================
3960 
3962 {
3963  // Generic four-particle correlation <exp[i(n1*phi1+n2*phi2+n3*phi3+n4*phi4)]>.
3964 
3965  TComplex four = Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n4,1)-Q(n1+n2,2)*Q(n3,1)*Q(n4,1)-Q(n2,1)*Q(n1+n3,2)*Q(n4,1)
3966  - Q(n1,1)*Q(n2+n3,2)*Q(n4,1)+2.*Q(n1+n2+n3,3)*Q(n4,1)-Q(n2,1)*Q(n3,1)*Q(n1+n4,2)
3967  + Q(n2+n3,2)*Q(n1+n4,2)-Q(n1,1)*Q(n3,1)*Q(n2+n4,2)+Q(n1+n3,2)*Q(n2+n4,2)
3968  + 2.*Q(n3,1)*Q(n1+n2+n4,3)-Q(n1,1)*Q(n2,1)*Q(n3+n4,2)+Q(n1+n2,2)*Q(n3+n4,2)
3969  + 2.*Q(n2,1)*Q(n1+n3+n4,3)+2.*Q(n1,1)*Q(n2+n3+n4,3)-6.*Q(n1+n2+n3+n4,4);
3970 
3971  return four;
3972 
3973 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::Four(Int_t n1, Int_t n2, Int_t n3, Int_t n4)
3974 
3975 //=======================================================================================================================
3976 
3978 {
3979  // Generic five-particle correlation <exp[i(n1*phi1+n2*phi2+n3*phi3+n4*phi4+n5*phi5)]>.
3980 
3981  TComplex five = Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n4,1)*Q(n5,1)-Q(n1+n2,2)*Q(n3,1)*Q(n4,1)*Q(n5,1)
3982  - Q(n2,1)*Q(n1+n3,2)*Q(n4,1)*Q(n5,1)-Q(n1,1)*Q(n2+n3,2)*Q(n4,1)*Q(n5,1)
3983  + 2.*Q(n1+n2+n3,3)*Q(n4,1)*Q(n5,1)-Q(n2,1)*Q(n3,1)*Q(n1+n4,2)*Q(n5,1)
3984  + Q(n2+n3,2)*Q(n1+n4,2)*Q(n5,1)-Q(n1,1)*Q(n3,1)*Q(n2+n4,2)*Q(n5,1)
3985  + Q(n1+n3,2)*Q(n2+n4,2)*Q(n5,1)+2.*Q(n3,1)*Q(n1+n2+n4,3)*Q(n5,1)
3986  - Q(n1,1)*Q(n2,1)*Q(n3+n4,2)*Q(n5,1)+Q(n1+n2,2)*Q(n3+n4,2)*Q(n5,1)
3987  + 2.*Q(n2,1)*Q(n1+n3+n4,3)*Q(n5,1)+2.*Q(n1,1)*Q(n2+n3+n4,3)*Q(n5,1)
3988  - 6.*Q(n1+n2+n3+n4,4)*Q(n5,1)-Q(n2,1)*Q(n3,1)*Q(n4,1)*Q(n1+n5,2)
3989  + Q(n2+n3,2)*Q(n4,1)*Q(n1+n5,2)+Q(n3,1)*Q(n2+n4,2)*Q(n1+n5,2)
3990  + Q(n2,1)*Q(n3+n4,2)*Q(n1+n5,2)-2.*Q(n2+n3+n4,3)*Q(n1+n5,2)
3991  - Q(n1,1)*Q(n3,1)*Q(n4,1)*Q(n2+n5,2)+Q(n1+n3,2)*Q(n4,1)*Q(n2+n5,2)
3992  + Q(n3,1)*Q(n1+n4,2)*Q(n2+n5,2)+Q(n1,1)*Q(n3+n4,2)*Q(n2+n5,2)
3993  - 2.*Q(n1+n3+n4,3)*Q(n2+n5,2)+2.*Q(n3,1)*Q(n4,1)*Q(n1+n2+n5,3)
3994  - 2.*Q(n3+n4,2)*Q(n1+n2+n5,3)-Q(n1,1)*Q(n2,1)*Q(n4,1)*Q(n3+n5,2)
3995  + Q(n1+n2,2)*Q(n4,1)*Q(n3+n5,2)+Q(n2,1)*Q(n1+n4,2)*Q(n3+n5,2)
3996  + Q(n1,1)*Q(n2+n4,2)*Q(n3+n5,2)-2.*Q(n1+n2+n4,3)*Q(n3+n5,2)
3997  + 2.*Q(n2,1)*Q(n4,1)*Q(n1+n3+n5,3)-2.*Q(n2+n4,2)*Q(n1+n3+n5,3)
3998  + 2.*Q(n1,1)*Q(n4,1)*Q(n2+n3+n5,3)-2.*Q(n1+n4,2)*Q(n2+n3+n5,3)
3999  - 6.*Q(n4,1)*Q(n1+n2+n3+n5,4)-Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n4+n5,2)
4000  + Q(n1+n2,2)*Q(n3,1)*Q(n4+n5,2)+Q(n2,1)*Q(n1+n3,2)*Q(n4+n5,2)
4001  + Q(n1,1)*Q(n2+n3,2)*Q(n4+n5,2)-2.*Q(n1+n2+n3,3)*Q(n4+n5,2)
4002  + 2.*Q(n2,1)*Q(n3,1)*Q(n1+n4+n5,3)-2.*Q(n2+n3,2)*Q(n1+n4+n5,3)
4003  + 2.*Q(n1,1)*Q(n3,1)*Q(n2+n4+n5,3)-2.*Q(n1+n3,2)*Q(n2+n4+n5,3)
4004  - 6.*Q(n3,1)*Q(n1+n2+n4+n5,4)+2.*Q(n1,1)*Q(n2,1)*Q(n3+n4+n5,3)
4005  - 2.*Q(n1+n2,2)*Q(n3+n4+n5,3)-6.*Q(n2,1)*Q(n1+n3+n4+n5,4)
4006  - 6.*Q(n1,1)*Q(n2+n3+n4+n5,4)+24.*Q(n1+n2+n3+n4+n5,5);
4007 
4008  return five;
4009 
4010 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::Five(Int_t n1, Int_t n2, Int_t n3, Int_t n4, Int_t n5)
4011 
4012 //=======================================================================================================================
4013 
4015 {
4016  // Generic six-particle correlation <exp[i(n1*phi1+n2*phi2+n3*phi3+n4*phi4+n5*phi5+n6*phi6)]>.
4017 
4018  TComplex six = Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n4,1)*Q(n5,1)*Q(n6,1)-Q(n1+n2,2)*Q(n3,1)*Q(n4,1)*Q(n5,1)*Q(n6,1)
4019  - Q(n2,1)*Q(n1+n3,2)*Q(n4,1)*Q(n5,1)*Q(n6,1)-Q(n1,1)*Q(n2+n3,2)*Q(n4,1)*Q(n5,1)*Q(n6,1)
4020  + 2.*Q(n1+n2+n3,3)*Q(n4,1)*Q(n5,1)*Q(n6,1)-Q(n2,1)*Q(n3,1)*Q(n1+n4,2)*Q(n5,1)*Q(n6,1)
4021  + Q(n2+n3,2)*Q(n1+n4,2)*Q(n5,1)*Q(n6,1)-Q(n1,1)*Q(n3,1)*Q(n2+n4,2)*Q(n5,1)*Q(n6,1)
4022  + Q(n1+n3,2)*Q(n2+n4,2)*Q(n5,1)*Q(n6,1)+2.*Q(n3,1)*Q(n1+n2+n4,3)*Q(n5,1)*Q(n6,1)
4023  - Q(n1,1)*Q(n2,1)*Q(n3+n4,2)*Q(n5,1)*Q(n6,1)+Q(n1+n2,2)*Q(n3+n4,2)*Q(n5,1)*Q(n6,1)
4024  + 2.*Q(n2,1)*Q(n1+n3+n4,3)*Q(n5,1)*Q(n6,1)+2.*Q(n1,1)*Q(n2+n3+n4,3)*Q(n5,1)*Q(n6,1)
4025  - 6.*Q(n1+n2+n3+n4,4)*Q(n5,1)*Q(n6,1)-Q(n2,1)*Q(n3,1)*Q(n4,1)*Q(n1+n5,2)*Q(n6,1)
4026  + Q(n2+n3,2)*Q(n4,1)*Q(n1+n5,2)*Q(n6,1)+Q(n3,1)*Q(n2+n4,2)*Q(n1+n5,2)*Q(n6,1)
4027  + Q(n2,1)*Q(n3+n4,2)*Q(n1+n5,2)*Q(n6,1)-2.*Q(n2+n3+n4,3)*Q(n1+n5,2)*Q(n6,1)
4028  - Q(n1,1)*Q(n3,1)*Q(n4,1)*Q(n2+n5,2)*Q(n6,1)+Q(n1+n3,2)*Q(n4,1)*Q(n2+n5,2)*Q(n6,1)
4029  + Q(n3,1)*Q(n1+n4,2)*Q(n2+n5,2)*Q(n6,1)+Q(n1,1)*Q(n3+n4,2)*Q(n2+n5,2)*Q(n6,1)
4030  - 2.*Q(n1+n3+n4,3)*Q(n2+n5,2)*Q(n6,1)+2.*Q(n3,1)*Q(n4,1)*Q(n1+n2+n5,3)*Q(n6,1)
4031  - 2.*Q(n3+n4,2)*Q(n1+n2+n5,3)*Q(n6,1)-Q(n1,1)*Q(n2,1)*Q(n4,1)*Q(n3+n5,2)*Q(n6,1)
4032  + Q(n1+n2,2)*Q(n4,1)*Q(n3+n5,2)*Q(n6,1)+Q(n2,1)*Q(n1+n4,2)*Q(n3+n5,2)*Q(n6,1)
4033  + Q(n1,1)*Q(n2+n4,2)*Q(n3+n5,2)*Q(n6,1)-2.*Q(n1+n2+n4,3)*Q(n3+n5,2)*Q(n6,1)
4034  + 2.*Q(n2,1)*Q(n4,1)*Q(n1+n3+n5,3)*Q(n6,1)-2.*Q(n2+n4,2)*Q(n1+n3+n5,3)*Q(n6,1)
4035  + 2.*Q(n1,1)*Q(n4,1)*Q(n2+n3+n5,3)*Q(n6,1)-2.*Q(n1+n4,2)*Q(n2+n3+n5,3)*Q(n6,1)
4036  - 6.*Q(n4,1)*Q(n1+n2+n3+n5,4)*Q(n6,1)-Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n4+n5,2)*Q(n6,1)
4037  + Q(n1+n2,2)*Q(n3,1)*Q(n4+n5,2)*Q(n6,1)+Q(n2,1)*Q(n1+n3,2)*Q(n4+n5,2)*Q(n6,1)
4038  + Q(n1,1)*Q(n2+n3,2)*Q(n4+n5,2)*Q(n6,1)-2.*Q(n1+n2+n3,3)*Q(n4+n5,2)*Q(n6,1)
4039  + 2.*Q(n2,1)*Q(n3,1)*Q(n1+n4+n5,3)*Q(n6,1)-2.*Q(n2+n3,2)*Q(n1+n4+n5,3)*Q(n6,1)
4040  + 2.*Q(n1,1)*Q(n3,1)*Q(n2+n4+n5,3)*Q(n6,1)-2.*Q(n1+n3,2)*Q(n2+n4+n5,3)*Q(n6,1)
4041  - 6.*Q(n3,1)*Q(n1+n2+n4+n5,4)*Q(n6,1)+2.*Q(n1,1)*Q(n2,1)*Q(n3+n4+n5,3)*Q(n6,1)
4042  - 2.*Q(n1+n2,2)*Q(n3+n4+n5,3)*Q(n6,1)-6.*Q(n2,1)*Q(n1+n3+n4+n5,4)*Q(n6,1)
4043  - 6.*Q(n1,1)*Q(n2+n3+n4+n5,4)*Q(n6,1)+24.*Q(n1+n2+n3+n4+n5,5)*Q(n6,1)
4044  - Q(n2,1)*Q(n3,1)*Q(n4,1)*Q(n5,1)*Q(n1+n6,2)+Q(n2+n3,2)*Q(n4,1)*Q(n5,1)*Q(n1+n6,2)
4045  + Q(n3,1)*Q(n2+n4,2)*Q(n5,1)*Q(n1+n6,2)+Q(n2,1)*Q(n3+n4,2)*Q(n5,1)*Q(n1+n6,2)
4046  - 2.*Q(n2+n3+n4,3)*Q(n5,1)*Q(n1+n6,2)+Q(n3,1)*Q(n4,1)*Q(n2+n5,2)*Q(n1+n6,2)
4047  - Q(n3+n4,2)*Q(n2+n5,2)*Q(n1+n6,2)+Q(n2,1)*Q(n4,1)*Q(n3+n5,2)*Q(n1+n6,2)
4048  - Q(n2+n4,2)*Q(n3+n5,2)*Q(n1+n6,2)-2.*Q(n4,1)*Q(n2+n3+n5,3)*Q(n1+n6,2)
4049  + Q(n2,1)*Q(n3,1)*Q(n4+n5,2)*Q(n1+n6,2)-Q(n2+n3,2)*Q(n4+n5,2)*Q(n1+n6,2)
4050  - 2.*Q(n3,1)*Q(n2+n4+n5,3)*Q(n1+n6,2)-2.*Q(n2,1)*Q(n3+n4+n5,3)*Q(n1+n6,2)
4051  + 6.*Q(n2+n3+n4+n5,4)*Q(n1+n6,2)-Q(n1,1)*Q(n3,1)*Q(n4,1)*Q(n5,1)*Q(n2+n6,2)
4052  + Q(n1+n3,2)*Q(n4,1)*Q(n5,1)*Q(n2+n6,2)+Q(n3,1)*Q(n1+n4,2)*Q(n5,1)*Q(n2+n6,2)
4053  + Q(n1,1)*Q(n3+n4,2)*Q(n5,1)*Q(n2+n6,2)-2.*Q(n1+n3+n4,3)*Q(n5,1)*Q(n2+n6,2)
4054  + Q(n3,1)*Q(n4,1)*Q(n1+n5,2)*Q(n2+n6,2)-Q(n3+n4,2)*Q(n1+n5,2)*Q(n2+n6,2)
4055  + Q(n1,1)*Q(n4,1)*Q(n3+n5,2)*Q(n2+n6,2)-Q(n1+n4,2)*Q(n3+n5,2)*Q(n2+n6,2)
4056  - 2.*Q(n4,1)*Q(n1+n3+n5,3)*Q(n2+n6,2)+Q(n1,1)*Q(n3,1)*Q(n4+n5,2)*Q(n2+n6,2)
4057  - Q(n1+n3,2)*Q(n4+n5,2)*Q(n2+n6,2)-2.*Q(n3,1)*Q(n1+n4+n5,3)*Q(n2+n6,2)
4058  - 2.*Q(n1,1)*Q(n3+n4+n5,3)*Q(n2+n6,2)+6.*Q(n1+n3+n4+n5,4)*Q(n2+n6,2)
4059  + 2.*Q(n3,1)*Q(n4,1)*Q(n5,1)*Q(n1+n2+n6,3)-2.*Q(n3+n4,2)*Q(n5,1)*Q(n1+n2+n6,3)
4060  - 2.*Q(n4,1)*Q(n3+n5,2)*Q(n1+n2+n6,3)-2.*Q(n3,1)*Q(n4+n5,2)*Q(n1+n2+n6,3)
4061  + 4.*Q(n3+n4+n5,3)*Q(n1+n2+n6,3)-Q(n1,1)*Q(n2,1)*Q(n4,1)*Q(n5,1)*Q(n3+n6,2)
4062  + Q(n1+n2,2)*Q(n4,1)*Q(n5,1)*Q(n3+n6,2)+Q(n2,1)*Q(n1+n4,2)*Q(n5,1)*Q(n3+n6,2)
4063  + Q(n1,1)*Q(n2+n4,2)*Q(n5,1)*Q(n3+n6,2)-2.*Q(n1+n2+n4,3)*Q(n5,1)*Q(n3+n6,2)
4064  + Q(n2,1)*Q(n4,1)*Q(n1+n5,2)*Q(n3+n6,2)-Q(n2+n4,2)*Q(n1+n5,2)*Q(n3+n6,2)
4065  + Q(n1,1)*Q(n4,1)*Q(n2+n5,2)*Q(n3+n6,2)-Q(n1+n4,2)*Q(n2+n5,2)*Q(n3+n6,2)
4066  - 2.*Q(n4,1)*Q(n1+n2+n5,3)*Q(n3+n6,2)+Q(n1,1)*Q(n2,1)*Q(n4+n5,2)*Q(n3+n6,2)
4067  - Q(n1+n2,2)*Q(n4+n5,2)*Q(n3+n6,2)-2.*Q(n2,1)*Q(n1+n4+n5,3)*Q(n3+n6,2)
4068  - 2.*Q(n1,1)*Q(n2+n4+n5,3)*Q(n3+n6,2)+6.*Q(n1+n2+n4+n5,4)*Q(n3+n6,2)
4069  + 2.*Q(n2,1)*Q(n4,1)*Q(n5,1)*Q(n1+n3+n6,3)-2.*Q(n2+n4,2)*Q(n5,1)*Q(n1+n3+n6,3)
4070  - 2.*Q(n4,1)*Q(n2+n5,2)*Q(n1+n3+n6,3)-2.*Q(n2,1)*Q(n4+n5,2)*Q(n1+n3+n6,3)
4071  + 4.*Q(n2+n4+n5,3)*Q(n1+n3+n6,3)+2.*Q(n1,1)*Q(n4,1)*Q(n5,1)*Q(n2+n3+n6,3)
4072  - 2.*Q(n1+n4,2)*Q(n5,1)*Q(n2+n3+n6,3)-2.*Q(n4,1)*Q(n1+n5,2)*Q(n2+n3+n6,3)
4073  - 2.*Q(n1,1)*Q(n4+n5,2)*Q(n2+n3+n6,3)+4.*Q(n1+n4+n5,3)*Q(n2+n3+n6,3)
4074  - 6.*Q(n4,1)*Q(n5,1)*Q(n1+n2+n3+n6,4)+6.*Q(n4+n5,2)*Q(n1+n2+n3+n6,4)
4075  - Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n5,1)*Q(n4+n6,2)+Q(n1+n2,2)*Q(n3,1)*Q(n5,1)*Q(n4+n6,2)
4076  + Q(n2,1)*Q(n1+n3,2)*Q(n5,1)*Q(n4+n6,2)+Q(n1,1)*Q(n2+n3,2)*Q(n5,1)*Q(n4+n6,2)
4077  - 2.*Q(n1+n2+n3,3)*Q(n5,1)*Q(n4+n6,2)+Q(n2,1)*Q(n3,1)*Q(n1+n5,2)*Q(n4+n6,2)
4078  - Q(n2+n3,2)*Q(n1+n5,2)*Q(n4+n6,2)+Q(n1,1)*Q(n3,1)*Q(n2+n5,2)*Q(n4+n6,2)
4079  - Q(n1+n3,2)*Q(n2+n5,2)*Q(n4+n6,2)-2.*Q(n3,1)*Q(n1+n2+n5,3)*Q(n4+n6,2)
4080  + Q(n1,1)*Q(n2,1)*Q(n3+n5,2)*Q(n4+n6,2)-Q(n1+n2,2)*Q(n3+n5,2)*Q(n4+n6,2)
4081  - 2.*Q(n2,1)*Q(n1+n3+n5,3)*Q(n4+n6,2)-2.*Q(n1,1)*Q(n2+n3+n5,3)*Q(n4+n6,2)
4082  + 6.*Q(n1+n2+n3+n5,4)*Q(n4+n6,2)+2.*Q(n2,1)*Q(n3,1)*Q(n5,1)*Q(n1+n4+n6,3)
4083  - 2.*Q(n2+n3,2)*Q(n5,1)*Q(n1+n4+n6,3)-2.*Q(n3,1)*Q(n2+n5,2)*Q(n1+n4+n6,3)
4084  - 2.*Q(n2,1)*Q(n3+n5,2)*Q(n1+n4+n6,3)+4.*Q(n2+n3+n5,3)*Q(n1+n4+n6,3)
4085  + 2.*Q(n1,1)*Q(n3,1)*Q(n5,1)*Q(n2+n4+n6,3)-2.*Q(n1+n3,2)*Q(n5,1)*Q(n2+n4+n6,3)
4086  - 2.*Q(n3,1)*Q(n1+n5,2)*Q(n2+n4+n6,3)-2.*Q(n1,1)*Q(n3+n5,2)*Q(n2+n4+n6,3)
4087  + 4.*Q(n1+n3+n5,3)*Q(n2+n4+n6,3)-6.*Q(n3,1)*Q(n5,1)*Q(n1+n2+n4+n6,4)
4088  + 6.*Q(n3+n5,2)*Q(n1+n2+n4+n6,4)+2.*Q(n1,1)*Q(n2,1)*Q(n5,1)*Q(n3+n4+n6,3)
4089  - 2.*Q(n1+n2,2)*Q(n5,1)*Q(n3+n4+n6,3)-2.*Q(n2,1)*Q(n1+n5,2)*Q(n3+n4+n6,3)
4090  - 2.*Q(n1,1)*Q(n2+n5,2)*Q(n3+n4+n6,3)+4.*Q(n1+n2+n5,3)*Q(n3+n4+n6,3)
4091  - 6.*Q(n2,1)*Q(n5,1)*Q(n1+n3+n4+n6,4)+6.*Q(n2+n5,2)*Q(n1+n3+n4+n6,4)
4092  - 6.*Q(n1,1)*Q(n5,1)*Q(n2+n3+n4+n6,4)+6.*Q(n1+n5,2)*Q(n2+n3+n4+n6,4)
4093  + 24.*Q(n5,1)*Q(n1+n2+n3+n4+n6,5)-Q(n1,1)*Q(n2,1)*Q(n3,1)*Q(n4,1)*Q(n5+n6,2)
4094  + Q(n1+n2,2)*Q(n3,1)*Q(n4,1)*Q(n5+n6,2)+Q(n2,1)*Q(n1+n3,2)*Q(n4,1)*Q(n5+n6,2)
4095  + Q(n1,1)*Q(n2+n3,2)*Q(n4,1)*Q(n5+n6,2)-2.*Q(n1+n2+n3,3)*Q(n4,1)*Q(n5+n6,2)
4096  + Q(n2,1)*Q(n3,1)*Q(n1+n4,2)*Q(n5+n6,2)-Q(n2+n3,2)*Q(n1+n4,2)*Q(n5+n6,2)
4097  + Q(n1,1)*Q(n3,1)*Q(n2+n4,2)*Q(n5+n6,2)-Q(n1+n3,2)*Q(n2+n4,2)*Q(n5+n6,2)
4098  - 2.*Q(n3,1)*Q(n1+n2+n4,3)*Q(n5+n6,2)+Q(n1,1)*Q(n2,1)*Q(n3+n4,2)*Q(n5+n6,2)
4099  - Q(n1+n2,2)*Q(n3+n4,2)*Q(n5+n6,2)-2.*Q(n2,1)*Q(n1+n3+n4,3)*Q(n5+n6,2)
4100  - 2.*Q(n1,1)*