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