AliPhysics  b0c77bb (b0c77bb)
AliAnalysisTaskTwoMultiCorrelations.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 // Analysis task for the data analysis of the correlations between the anisotropic flow //
18 // harmonics v_n with the Pb-Pb data taken by the ALICE experiment. //
19 // The current script computes the multiparticle correlators using the method of the //
20 // Q-vectors for a maximum of 6 different harmonics and 8 particles). //
21 // //
22 // Author: Cindy Mordasini (cindy.mordasini@cern.ch) //
23 //--------------------------------------------------------------------------------------//
24 
25 #include "Riostream.h"
26 #include <vector>
28 #include "AliLog.h"
29 #include "AliAODEvent.h"
30 #include "AliAODInputHandler.h"
31 #include "AliMCEvent.h"
32 #include "AliMCEventHandler.h"
33 #include "AliAnalysisManager.h"
34 #include "AliMultSelection.h"
35 #include "TList.h"
36 #include "TH1D.h"
37 #include "TFile.h"
38 #include "TComplex.h"
39 #include "TMath.h"
40 
41 using std::cout;
42 using std::endl;
43 
45 
46 //======================================================================================//
47 
49  AliAnalysisTaskSE(name),
50  fMaxNumberCorrelations(8),
51  fMaxFlowHarmonic(6),
52  fNumberHarmonicsInSC(4),
53  fAnalysisType(NULL),
54  fProcessBothKineAndReco(kFALSE),
55  fProcessOnlyKine(kFALSE),
56  fProcessOnlyReco(kFALSE),
57  fUseParticleWeights(kFALSE),
58  fComputeNestedLoops(kFALSE),
59  fMainList(NULL),
60  fControlListPreCuts(NULL),
61  fControlListPostCuts(NULL),
62  fFinalList(NULL),
63  fListTwoParticles(NULL),
64  fListThreeParticles(NULL),
65  fListFourParticles(NULL),
66  fListSixParticles(NULL),
67  fCentralityMin(0.),
68  fCentralityMax(100.),
69  fPtMin(0.2),
70  fPtMax(5.),
71  fEtaMin(-0.8),
72  fEtaMax(0.8),
73  fHarmonicOne(2),
74  fHarmonicTwo(-2),
75  fHarmonicThree(3),
76  fHarmonicFour(-3),
77  fHarmonicFive(4),
78  fHarmonicSix(-4),
79  fHarmonicSeven(5),
80  fHarmonicEight(-5),
81  fHistoCentralityPreCuts(NULL),
82  fHistoPtPreCuts(NULL),
83  fHistoEtaPreCuts(NULL),
84  fHistoPhiPreCuts(NULL),
85  fHistoMultiplicityPreCuts(NULL),
86  fHistoMultiplicityPostCuts(NULL),
87  fHistoPtPostCuts(NULL),
88  fHistoEtaPostCuts(NULL),
89  fHistoPhiPostCuts(NULL),
90  fProfileCosineEightParticles(NULL)
91 {
92 // Constructor of the class.
93  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
94 
95 // Create the main list.
96  fMainList = new TList();
97  fMainList->SetName("outputAnalysis");
98  fMainList->SetOwner(kTRUE);
99 
100 // Define the input and output slots.
101  DefineOutput(1, TList::Class());
102 
103 // Initialise the fQvectors to zero.
104  InitialiseArraysOfQvectors();
105 
106 // Initialise the pointers of the TProfiles to NULL.
107  InitialiseArraysOfTProfiles();
108 
109 // Use the particle weights?
110  if(useParticleWeights)
111  {
112  // Not needed for 2010 Heavy Ions data.
113  // TBI for later data taking periods...
114  } // End: if(useParticleWeights)
115 } // End: AliAnalysisTaskTwoMultiCorrelations(const chat*, Bool_t)
116 
117 //======================================================================================//
118 
121  fMaxNumberCorrelations(8),
122  fMaxFlowHarmonic(6),
123  fNumberHarmonicsInSC(4),
124  fAnalysisType(NULL),
125  fProcessBothKineAndReco(kFALSE),
126  fProcessOnlyKine(kFALSE),
127  fProcessOnlyReco(kFALSE),
128  fUseParticleWeights(kFALSE),
129  fComputeNestedLoops(kFALSE),
130  fMainList(NULL),
131  fControlListPreCuts(NULL),
132  fControlListPostCuts(NULL),
133  fFinalList(NULL),
134  fListTwoParticles(NULL),
135  fListThreeParticles(NULL),
136  fListFourParticles(NULL),
137  fListSixParticles(NULL),
138  fCentralityMin(0.),
139  fCentralityMax(100.),
140  fPtMin(0.2),
141  fPtMax(5.),
142  fEtaMin(-0.8),
143  fEtaMax(0.8),
144  fHarmonicOne(2),
145  fHarmonicTwo(-2),
146  fHarmonicThree(3),
147  fHarmonicFour(-3),
148  fHarmonicFive(4),
149  fHarmonicSix(-4),
150  fHarmonicSeven(5),
151  fHarmonicEight(-5),
152  fHistoCentralityPreCuts(NULL),
153  fHistoPtPreCuts(NULL),
154  fHistoEtaPreCuts(NULL),
155  fHistoPhiPreCuts(NULL),
156  fHistoMultiplicityPreCuts(NULL),
157  fHistoMultiplicityPostCuts(NULL),
158  fHistoPtPostCuts(NULL),
159  fHistoEtaPostCuts(NULL),
160  fHistoPhiPostCuts(NULL),
161  fProfileCosineEightParticles(NULL)
162 {
163 // Dummy constructor of the class.
164  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
165 
166 // Initialise the fQvectors to zero.
168 
169 // Initialise the pointers of the TProfiles to NULL.
171 } // End: AliAnalysisTaskTwoMultiCorrelations()
172 
173 //======================================================================================//
174 
176 {
177 // Destructor of the class.
178  if(fMainList) {delete fMainList;}
179 } // End: ~AliAnalysisTaskTwoMultiCorrelations()
180 
181 //======================================================================================//
182 
184 {
185 // Initialisation of the lists at the beginning of the analysis.
186 // Avoid name clashes.
187  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
188  TH1::AddDirectory(kFALSE);
189 
190 // Book and nest all the lists.
191  this->BookAndNestAllLists();
192 
193 // Book all the secondary lists.
194  this->BookControlListPreCuts();
195  this->BookControlListPostCuts();
196  this->BookFinalList();
197 
198 // Initialise the pointer for the analysis type.
199  fAnalysisType = new TString();
200 
201 // Still avoid name clashes.
202  TH1::AddDirectory(oldHistAddStatus);
203  PostData(1,fMainList);
204 } // End: UserCreateOutputObjects()
205 
206 //======================================================================================//
207 
209 {
210 // Do the main analysis for each event.
211 // TBI: find a way to include non-unit particle weights and select them via fUseParticleWeights...
212  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::UserExec(Option_t *)";
213 
214 // Select of the analysis type (MC/AOD) from TaskSE.
215  AliMCEvent *currentMCEvent = MCEvent();
216  AliAODEvent *currentAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
217  //if(!currentAODEvent){return;} // Protection against NULL pointer
219  {
220  Fatal(sMethodName.Data(),"ERROR: only one fProcess must be kTRUE in 'SetAnalysisType'");
221  }
222  else if (fProcessBothKineAndReco) {*fAnalysisType="MC_AOD";}
223  else if (fProcessOnlyKine) {*fAnalysisType="MC";}
224  else if (fProcessOnlyReco) {*fAnalysisType="AOD";}
225 
226 // Call of the specific analysis.
227  if (fAnalysisType->EqualTo("AOD")) {AODanalysis(currentAODEvent);}
228  else if (fAnalysisType->EqualTo("MC")) {MCanalysis(currentMCEvent);}
229  else if (fAnalysisType->EqualTo("MC_AOD")) {Fatal(sMethodName.Data(),"TBI...");}
230  else {Fatal(sMethodName.Data(),"ERROR: fAnalysisType not defined");}
231 
232 // ...
233 
234 // PostData.
235  PostData(1, fMainList);
236 } // End: UserExec(Option_t *)
237 
238 //======================================================================================//
239 
241 {
242 // Post-analysis done once all the events are done.
243 // Access the merged main list.
244  fMainList = (TList*)GetOutputData(1);
245  if(!fMainList){exit(1);}
246 
247 // Create the output file and save the main list in it
248  TFile *outputFile = new TFile("AnalysisResults.root", "RECREATE");
249  fMainList->Write(fMainList->GetName(),TObject::kSingleKey);
250  delete outputFile;
251 } // End: Terminate(Option_t *)
252 
253 //======================================================================================//
254 // Methods called in the constructor.
255 //--------------------------------------------------------------------------------------//
256 
258 {
259 // Initialisation of the Q-vectors to zero.
260  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
261  {
262  for (Int_t iPower = 0; iPower < 9; iPower++)
263  {
264  fQvectors[iHarmo][iPower] = TComplex(0.,0.);
265  } // End: for (Int_t iPower = 0; iPower < 9; iPower++)
266  } // End: for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
267 } // End: InitialiseArraysOfQvectors()
268 
269 //--------------------------------------------------------------------------------------//
270 
272 {
273 // Initialisation of the pointers to the TProfiles.
274  for (Int_t i = 0; i < 6; i++)
275  {
276  fProfileCosineTwoParticles[i] = NULL;
278  fProfileCosineFourParticles[i] = NULL;
280  } // End: for (Int_t i = 0; i < 6; i++)
281 
282  for (Int_t j = 0; j < 5; j++)
283  {
284  fProfileTwoCosine[j] = NULL;
286  } // End: for (Int_t j = 0; j < 5; j++)
287 
288  for (Int_t k = 0; k < 4; k++)
289  {
292  fProfileCosineSixParticles[k] = NULL;
293  } // End: for (Int_t k = 0; k < 4; k++)
294 
295 } // End: InitialiseArraysOfTProfiles()
296 
297 //======================================================================================//
298 // Methods called in UserCreateOutputObjects.
299 //--------------------------------------------------------------------------------------//
300 
302 {
303 // Book all the lists.
304  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::BookAndNestAllLists()";
305  if(!fMainList){Fatal(sMethodName.Data(),"Main list 'fMainList' is NULL");}
306 
307 // Control histograms before the cuts.
308  fControlListPreCuts = new TList();
309  fControlListPreCuts->SetName("ControlListPreCuts");
310  fControlListPreCuts->SetOwner(kTRUE);
311  fMainList->Add(fControlListPreCuts);
312 
313 // Control histograms after the cuts
314  fControlListPostCuts = new TList();
315  fControlListPostCuts->SetName("ControlListPostCuts");
316  fControlListPostCuts->SetOwner(kTRUE);
318 
319 // Results for the multiparticle correlations divided by number of particles involved.
320  fFinalList = new TList();
321  fFinalList->SetName("FinalList");
322  fFinalList->SetOwner(kTRUE);
323  fMainList->Add(fFinalList);
324 
325  fListTwoParticles = new TList();
326  fListTwoParticles->SetName("ListTwoParticles");
327  fListTwoParticles->SetOwner(kTRUE);
329 
330  fListThreeParticles = new TList();
331  fListThreeParticles->SetName("ListThreeParticles");
332  fListThreeParticles->SetOwner(kTRUE);
334 
335  fListFourParticles = new TList();
336  fListFourParticles->SetName("ListFourParticles");
337  fListFourParticles->SetOwner(kTRUE);
339 
340  fListSixParticles = new TList();
341  fListSixParticles->SetName("ListSixParticles");
342  fListSixParticles->SetOwner(kTRUE);
344 } // End: BookAndNestAllLists()
345 
346 //--------------------------------------------------------------------------------------//
347 
349 {
350 // Book the histograms with the observables before the cuts.
351 // Centrality distribution.
352  fHistoCentralityPreCuts = new TH1D("fHistoCentralityPreCuts", "Centrality distribution before the cuts", 100, 0., 100.);
353  fHistoCentralityPreCuts->SetStats(kTRUE);
354  fHistoCentralityPreCuts->GetXaxis()->SetTitle("Centrality percentile");
356 
357 // Transverse momentum distribution.
358  fHistoPtPreCuts = new TH1D("fHistoPtPreCuts", "Transverse momentum distribution before the cuts", 1000, 0., 20.);
359  fHistoPtPreCuts->SetStats(kTRUE);
360  fHistoPtPreCuts->GetXaxis()->SetTitle("p_{T}");
362 
363 // Pseudorapidity distribution.
364  fHistoEtaPreCuts = new TH1D("fHistoEtaPreCuts", "Pseudorapidity distribution before the cuts", 1000, -1., 1.);
365  fHistoEtaPreCuts->SetStats(kTRUE);
366  fHistoEtaPreCuts->GetXaxis()->SetTitle("#eta");
368 
369 // Azimuthal angles distribution.
370  fHistoPhiPreCuts = new TH1D("fHistoPhiPreCuts", "Azimuthal angles distribution before the cuts", 1000, 0., 6.3);
371  fHistoPhiPreCuts->SetStats(kTRUE);
372  fHistoPhiPreCuts->GetXaxis()->SetTitle("#phi");
374 
375 // Multiplicity distribution.
376  fHistoMultiplicityPreCuts = new TH1D("fHistoMultiplicityPreCuts", "Multiplicity distribution before the cuts", 10000000, 0., 10000000.);
377  fHistoMultiplicityPreCuts->SetStats(kTRUE);
378  fHistoMultiplicityPreCuts->GetXaxis()->SetTitle("Multiplicity");
380 } // End: BookControlListPreCuts()
381 
382 //--------------------------------------------------------------------------------------//
383 
385 {
386 // Book the histograms with the observables after the cuts.
387 // Multiplicity distribution.
388  fHistoMultiplicityPostCuts = new TH1D("fHistoMultiplicityPostCuts", "Multiplicity distribution after the cuts", 10000000, 0., 10000000.);
389  fHistoMultiplicityPostCuts->SetStats(kTRUE);
390  fHistoMultiplicityPostCuts->GetXaxis()->SetTitle("Multiplicity");
392 
393 // Transverse momentum distribution.
394  fHistoPtPostCuts = new TH1D("fHistoPtPostCuts", "Transverse momentum distribution before the cuts", 1000, 0., 20.);
395  fHistoPtPostCuts->SetStats(kTRUE);
396  fHistoPtPostCuts->GetXaxis()->SetTitle("p_{T}");
398 
399 // Pseudorapidity distribution.
400  fHistoEtaPostCuts = new TH1D("fHistoEtaPostCuts", "Pseudorapidity distribution before the cuts", 1000, -1., 1.);
401  fHistoEtaPostCuts->SetStats(kTRUE);
402  fHistoEtaPostCuts->GetXaxis()->SetTitle("#eta");
404 
405 // Azimuthal angles distribution.
406  fHistoPhiPostCuts = new TH1D("fHistoPhiPostCuts", "Azimuthal angles distribution before the cuts", 1000, 0., 6.3);
407  fHistoPhiPostCuts->SetStats(kTRUE);
408  fHistoPhiPostCuts->GetXaxis()->SetTitle("#phi");
410 } // End: BookControlListPostCuts()
411 
412 //--------------------------------------------------------------------------------------//
413 
415 {
416 // Book the TProfiles with the results for the multiparticle correlations.
417 // For 2-particle correlations.
418  for (Int_t t = 0; t < 6; t++)
419  {
420  fProfileCosineTwoParticles[t] = new TProfile("", "", 1, 0., 1.);
421  fProfileCosineTwoParticles[t]->SetStats(kTRUE);
422  fProfileCosineTwoParticles[t]->Sumw2();
423  fProfileCosineTwoParticles[t]->SetName(Form("fProfileCosineTwoParticles_%d", t));
425 
427  {
428  fProfileCosineTwoNestedLoops[t] = new TProfile("", "", 1, 0., 1.);
429  fProfileCosineTwoNestedLoops[t]->SetStats(kTRUE);
430  fProfileCosineTwoNestedLoops[t]->Sumw2();
431  fProfileCosineTwoNestedLoops[t]->SetName(Form("fProfileCosineTwoNestedLoops_%d", t));
433  } // End: if (fComputeNestedLoops)
434  } // End: for (Int_t t = 0; t < 6; t++)
435 
436  for (Int_t c = 0; c < 5; c++)
437  {
438  fProfileTwoCosine[c] = new TProfile("", "", 1, 0., 1.);
439  fProfileTwoCosine[c]->SetStats(kTRUE);
440  fProfileTwoCosine[c]->Sumw2();
441  fProfileTwoCosine[c]->SetName(Form("fProfileTwoCosine_%d", c));
443 
445  {
446  fProfileTwoCosineNestedLoops[c] = new TProfile("", "", 1, 0., 1.);
447  fProfileTwoCosineNestedLoops[c]->SetStats(kTRUE);
449  fProfileTwoCosineNestedLoops[c]->SetName(Form("fProfileTwoCosineNestedLoops_%d", c));
451  } // End: if (fComputeNestedLoops)
452  } // End: for (Int_t c = 0; c < 4; c++)
453 
454 // For 3-particle correlations.
455  for (Int_t th = 0; th < 4; th++)
456  {
457  fProfileCosineThreeParticles[th] = new TProfile("","", 1, 0., 1.);
458  fProfileCosineThreeParticles[th]->SetStats(kTRUE);
459  fProfileCosineThreeParticles[th]->Sumw2();
460  fProfileCosineThreeParticles[th]->SetName(Form("fProfileCosineThreeParticles_%d", th));
462 
464  {
465  fProfileCosineThreeNestedLoops[th] = new TProfile("","", 1, 0., 1.);
466  fProfileCosineThreeNestedLoops[th]->SetStats(kTRUE);
467  fProfileCosineThreeNestedLoops[th]->Sumw2();
468  fProfileCosineThreeNestedLoops[th]->SetName(Form("fProfileCosineThreeNestedLoops_%d", th));
470  } // End: if (fComputeNestedLoops)
471  } // End: for (Int_t th = 0; th < 4; th++)
472 
473 // For 4-particle correlations.
474  for (Int_t f = 0; f < 6; f++)
475  {
476  fProfileCosineFourParticles[f] = new TProfile("", "", 1, 0., 1.);
477  fProfileCosineFourParticles[f]->SetStats(kTRUE);
478  fProfileCosineFourParticles[f]->Sumw2();
479  fProfileCosineFourParticles[f]->SetName(Form("fProfileCosineFourParticles_%d", f));
481 
483  {
484  fProfileCosineFourNestedLoops[f] = new TProfile("", "", 1, 0., 1.);
485  fProfileCosineFourNestedLoops[f]->SetStats(kTRUE);
486  fProfileCosineFourNestedLoops[f]->Sumw2();
487  fProfileCosineFourNestedLoops[f]->SetName(Form("fProfileCosineFourNestedLoops_%d", f));
489  } // End: if (fComputeNestedLoops)
490  } // End: for (Int_t f = 0; f < 6; f++)
491 
492 // For 6- and 8-particle correlations.
493  for (Int_t s = 0; s < 4; s++)
494  {
495  fProfileCosineSixParticles[s] = new TProfile("", "", 1, 0., 1.);
496  fProfileCosineSixParticles[s]->SetStats(kTRUE);
497  fProfileCosineSixParticles[s]->Sumw2();
498  fProfileCosineSixParticles[s]->SetName(Form("fProfileCosineSixParticles_%d", s));
500  } // End: for (Int_t s = 0; s < 4; s++)
501 
502  fProfileCosineEightParticles = new TProfile("", "", 1, 0., 1.);
503  fProfileCosineEightParticles->SetStats(kTRUE);
505  fProfileCosineEightParticles->SetName("fProfileCosineEightParticles");
507 
508 } // End: BookFinalList()
509 
510 //======================================================================================//
511 // Methods called in UserExec(Option_t *).
512 //--------------------------------------------------------------------------------------//
513 
515 {
516 // Compute the multi-particle correlations for the given AOD event.
517 // Define the centrality according to ALICE standards.
518  AliMultSelection *ams = (AliMultSelection*)aAODevent->FindListObject("MultSelection");
519  if(!ams){return;}
520  if(ams->GetMultiplicityPercentile("V0M") >= fCentralityMin && ams->GetMultiplicityPercentile("V0M") < fCentralityMax)
521  {
522  fHistoCentralityPreCuts->Fill(ams->GetMultiplicityPercentile("V0M"));
523  }
524  else {return;} // Event not belong to the specified centrality class.
525 
526 //**************************************************************************************//
527 // Determine which particles pass the cuts.
528  long long nParticlesBeforeCuts = aAODevent->GetNumberOfTracks(); // Number of particles before the cuts.
529  long long nParticlesAfterCuts = 0; // Number of particles after the cuts.
530  Int_t *particlesIndices = new Int_t[nParticlesBeforeCuts](); // List of the particles saved with the corresponding index = 1 if they pass the cuts or 0 if they do not.
531 
532  for (long long iParticle = 0; iParticle < nParticlesBeforeCuts; iParticle++)
533  {
534  // Pointer to the current particle.
535  AliAODTrack *currentParticle = dynamic_cast<AliAODTrack*>(aAODevent->GetTrack(iParticle));
536  if(!currentParticle){continue;} // Protection against NULL pointer.
537  if(!currentParticle->TestFilterBit(128)){continue;} // Filter bit 128 denotes TPC-only tracks.
538 
539  // Control histograms of some observables before the cuts.
540  Double_t ptPreCuts = currentParticle->Pt(); // Transverse momentum.
541  fHistoPtPreCuts->Fill(ptPreCuts);
542  Double_t etaPreCuts = currentParticle->Eta(); // Pseudorapidity.
543  fHistoEtaPreCuts->Fill(etaPreCuts);
544  Double_t phiPreCuts = currentParticle->Phi(); // Azimuthal angles.
545  fHistoPhiPreCuts->Fill(phiPreCuts);
546 
547  // Apply the cuts and set the flag 1/0 of particlesIndices[].
548  if ( (fEtaMin < etaPreCuts) && (etaPreCuts < fEtaMax) && (fPtMin < ptPreCuts) && (ptPreCuts < fPtMax) )
549  {
550  nParticlesAfterCuts++;
551  particlesIndices[iParticle] = 1;
552  } // End: if ("cuts")
553  else {particlesIndices[iParticle] = 0;}
554  } // End: for (long long iParticle = 0; iParticle < nParticlesBeforeCuts; iParticle++)
555 
556 //**************************************************************************************//
557 // Control histograms for the multiplicity before and after the cuts.
558  fHistoMultiplicityPreCuts->Fill(nParticlesBeforeCuts);
559  fHistoMultiplicityPostCuts->Fill(nParticlesAfterCuts);
560 
561 //**************************************************************************************//
562 // Prepare the variables for the AOD analysis.
563  Double_t *pt = new Double_t[nParticlesAfterCuts](); // Transverse momenta.
564  Double_t *eta = new Double_t[nParticlesAfterCuts](); // Pseudorapidity.
565  Double_t *phi = new Double_t[nParticlesAfterCuts](); // Azimuthal angles.
566  Double_t *particleWeights = new Double_t[nParticlesAfterCuts](); // Particle weights
567  Int_t index = 0; // Index of the particles which pass the cuts.
568 
569  if (nParticlesAfterCuts <= 2*fNumberHarmonicsInSC) {return;} // Analysis done only if the number of particles after the cuts is larger than the maximum number of particles needed for the studied GSC.
570  for (long long iiParticle = 0; iiParticle < nParticlesBeforeCuts; iiParticle++)
571  {
572  // Pointer to the current particle.
573  AliAODTrack *keptParticle = dynamic_cast<AliAODTrack*>(aAODevent->GetTrack(iiParticle));
574  if(!keptParticle){continue;} // Protection against NULL pointer.
575  if(!keptParticle->TestFilterBit(128)){continue;} // Filter bit 128 denotes TPC-only tracks.
576 
577  // Control histograms.
578  if (particlesIndices[iiParticle] == 1)
579  {
580  pt[index] = keptParticle->Pt(); // Transverse momentum.
581  fHistoPtPostCuts->Fill(pt[index]);
582  eta[index] = keptParticle->Eta();
583  fHistoEtaPostCuts->Fill(eta[index]);
584  phi[index] = keptParticle->Phi();
585  fHistoPhiPostCuts->Fill(phi[index]);
586  particleWeights[index] = 1.;
587  //if (fUseParticleWeights) {continue;} // TBI
588  //else {particleWeights[iiParticle] = 1.;}
589  index++;
590  } // End: if (particlesIndices[iiParticle] == 1)
591  else {continue;}
592  } // End: for (long long iiParticle = 0; iiParticle < nParticlesBeforeCuts; iiParticle++)
593 
594 //**************************************************************************************//
595 // Calculate all the Q-vectors for this event.
596  CalculateQvectors(nParticlesAfterCuts, phi, particleWeights);
597 
598 //**************************************************************************************//
599 // Do the full analysis for this event.
600  GSCfullAnalysis(nParticlesAfterCuts, phi, particleWeights);
601 
602 //**************************************************************************************//
603 // Resets of the general variables of the events.
604  delete [] pt;
605  delete [] eta;
606  delete [] phi;
607  delete [] particleWeights;
608  index = 0;
609 } // End: AODanalysis(AliAODEvent*)
610 
611 //--------------------------------------------------------------------------------------//
613 {
614 // Compute the multi-particle correlations for the given MC event.
615 // Define the centrality according to ALICE standards in the case of reconstructed particles.
616  if (!fProcessOnlyKine)
617  {
618  AliMultSelection *ams = (AliMultSelection*)aMCevent->FindListObject("MultSelection");
619  if(!ams){return;}
620  if(ams->GetMultiplicityPercentile("V0M") >= fCentralityMin && ams->GetMultiplicityPercentile("V0M") < fCentralityMax)
621  {
622  fHistoCentralityPreCuts->Fill(ams->GetMultiplicityPercentile("V0M"));
623  }
624  else {return;} // Event not belong to the specified centrality class.
625  } // End: if (!fProcessOnlyKine)
626 
627 // Determine which particles pass the cuts.
628  long long nParticlesBeforeCuts = aMCevent->GetNumberOfTracks(); // Number of particles before the cuts.
629  long long nParticlesAfterCuts = 0; // Number of particles after the cuts.
630  Int_t *particlesIndices = new Int_t[nParticlesBeforeCuts](); // List of the particles saved with the corresponding index = 1 if they pass the cuts or 0 if they do not.
631 
632  for (long long iParticle = 0; iParticle < nParticlesBeforeCuts; iParticle++)
633  {
634  // Pointer to the current particle.
635  AliAODMCParticle *currentParticle = dynamic_cast<AliAODMCParticle*>(aMCevent->GetTrack(iParticle));
636  if(!currentParticle){continue;} // Protection against NULL pointer.
637 
638  // Control histograms of some observables before the cuts.
639  Double_t ptPreCuts = currentParticle->Pt(); // Transverse momentum.
640  fHistoPtPreCuts->Fill(ptPreCuts);
641  Double_t etaPreCuts = currentParticle->Eta(); // Pseudorapidity.
642  fHistoEtaPreCuts->Fill(etaPreCuts);
643  Double_t phiPreCuts = currentParticle->Phi(); // Azimuthal angles.
644  fHistoPhiPreCuts->Fill(phiPreCuts);
645 
646  // Apply the cuts and set the flag 1/0 of particlesIndices[].
647  if ( (fEtaMin < etaPreCuts) && (etaPreCuts < fEtaMax) && (fPtMin < ptPreCuts) && (ptPreCuts < fPtMax) )
648  {
649  nParticlesAfterCuts++;
650  particlesIndices[iParticle] = 1;
651  } // End: if ("cuts")
652  else {particlesIndices[iParticle] = 0;}
653  } // End: for (long long iParticle = 0; iParticle < nParticlesBeforeCuts; iParticle++)
654 
655 //**************************************************************************************//
656 // Control histograms for the multiplicity before and after the cuts.
657  fHistoMultiplicityPreCuts->Fill(nParticlesBeforeCuts);
658  fHistoMultiplicityPostCuts->Fill(nParticlesAfterCuts);
659 
660 //**************************************************************************************//
661 // Prepare the variables for the MC analysis.
662  Double_t *pt = new Double_t[nParticlesAfterCuts](); // Transverse momenta.
663  Double_t *eta = new Double_t[nParticlesAfterCuts](); // Pseudorapidity.
664  Double_t *phi = new Double_t[nParticlesAfterCuts](); // Azimuthal angles.
665  Double_t *particleWeights = new Double_t[nParticlesAfterCuts](); // Particle weights
666  Int_t index = 0; // Index of the particles which pass the cuts.
667 
668  if (nParticlesAfterCuts <= 2*fNumberHarmonicsInSC) {return;} // Analysis done only if the number of particles after the cuts is larger than the maximum number of particles needed for the studied GSC.
669  for (long long iiParticle = 0; iiParticle < nParticlesBeforeCuts; iiParticle++)
670  {
671  // Pointer to the current particle.
672  AliAODMCParticle *keptParticle = dynamic_cast<AliAODMCParticle*>(aMCevent->GetTrack(iiParticle));
673  if(!keptParticle){continue;} // Protection against NULL pointer.
674 
675  // Control histograms.
676  if (particlesIndices[iiParticle] == 1)
677  {
678  pt[index] = keptParticle->Pt(); // Transverse momentum.
679  fHistoPtPostCuts->Fill(pt[index]);
680  eta[index] = keptParticle->Eta();
681  fHistoEtaPostCuts->Fill(eta[index]);
682  phi[index] = keptParticle->Phi();
683  fHistoPhiPostCuts->Fill(phi[index]);
684  particleWeights[index] = 1.;
685  //if (fUseParticleWeights) {continue;} // TBI
686  //else {particleWeights[iiParticle] = 1.;}
687  index++;
688  } // End: if (particlesIndices[iiParticle] == 1)
689  else {continue;}
690  } // End: for (long long iiParticle = 0; iiParticle < nParticlesBeforeCuts; iiParticle++)
691 
692 //**************************************************************************************//
693 // Calculate all the Q-vectors for this event.
694  CalculateQvectors(nParticlesAfterCuts, phi, particleWeights);
695 
696 //**************************************************************************************//
697 // Do the full analysis for this event.
698  GSCfullAnalysis(nParticlesAfterCuts, phi, particleWeights);
699 
700 //**************************************************************************************//
701 // Resets of the general variables of the events.
702  delete [] pt;
703  delete [] eta;
704  delete [] phi;
705  delete [] particleWeights;
706  index = 0;
707 } // End: MCDanalysis(AliMCEvent*)
708 //--------------------------------------------------------------------------------------//
709 
710 void AliAnalysisTaskTwoMultiCorrelations::CalculateQvectors(long long nParticles, Double_t angles[], Double_t weights[])
711 {
712 // Calculate all the Q-vectors for a given set of azimuthal angles and particles weights.
713  Double_t iAngle = 0.; // Angle of the current particle.
714  Double_t iWeight = 1.; // Particle weight of the current particle.
715  Double_t weightToPowerP = 1.; // Particle weight rised to the power p.
716 
717 // Ensure all the Q-vectors are initially zero.
718  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
719  {
720  for (Int_t iPower = 0; iPower < 9; iPower++)
721  {
722  fQvectors[iHarmo][iPower] = TComplex(0.,0.);
723  } // End: for (Int_t iPower = 0; iPower < 9; iPower++)
724  } // End: for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
725 
726 // Compute the Q-vectors.
727  for (long long iParticle = 0; iParticle < nParticles; iParticle++)
728  {
729  iAngle = angles[iParticle];
730  if (fUseParticleWeights) {iWeight = weights[iParticle];}
731  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
732  {
733  for (Int_t iPower = 0; iPower < 9; iPower++)
734  {
735  if (fUseParticleWeights) {weightToPowerP = TMath::Power(iWeight, iPower);}
736  fQvectors[iHarmo][iPower] += TComplex(weightToPowerP*TMath::Cos(iHarmo*iAngle), weightToPowerP*TMath::Sin(iHarmo*iAngle));
737  } // End: for (Int_t iPower = 0; iPower < 9; iPower++)
738  } // End: for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
739  } // End: for (Int_t iParticle = 0; iParticle < nParticles; iParticle++)
740 } // End: CalculateQvectors(long long, Double_t[], Double_t[])
741 
742 //--------------------------------------------------------------------------------------//
743 
744 void AliAnalysisTaskTwoMultiCorrelations::GSCfullAnalysis(long long nParticles, Double_t angles[], Double_t weights[])
745 {
746 // Compute all the needed multiparticle correlators (with possible cross-check from nested loops) for the chosen harmonics.
747  Int_t harmonicsArray[8] = {fHarmonicOne, fHarmonicTwo, fHarmonicThree, fHarmonicFour, fHarmonicFive, fHarmonicSix, fHarmonicSeven, fHarmonicEight}; // Harmonics (n_1,... n_8).
748 
749 // Compute the denominator of each case of correlators.
750  Int_t twoZerosArray[2] = {0}; // For 2p correlations.
751  Double_t twoParticleDenominator = CalculateRecursion(2, twoZerosArray).Re();
752  Int_t threeZerosArray[3] = {0}; // For 3p correlations.
753  Double_t threeParticleDenominator = CalculateRecursion(3, threeZerosArray).Re();
754  Int_t fourZerosArray[4] = {0}; // For 4p correlations.
755  Double_t fourParticleDenominator = CalculateRecursion(4, fourZerosArray).Re();
756  Int_t sixZerosArray[6] = {0}; // For 6p correlations.
757  Double_t sixParticleDenominator = CalculateRecursion(6, sixZerosArray).Re();
758  Int_t eightZerosArray[8] = {0}; // For 8p correlations.
759  Double_t eightParticleDenominator = CalculateRecursion(8, eightZerosArray).Re();
760  eightParticleDenominator++; // Dummy line since this variable is not useful yet.
761 
762 //**************************************************************************************//
763 // Compute the 2-particle correlations.
764  Int_t firstHarmonicArray[2] = {harmonicsArray[0], harmonicsArray[1]}; // Array with (k,-k).
765  Int_t secondHarmonicArray[2] = {harmonicsArray[2], harmonicsArray[3]}; // Array with (l,-l).
766  Int_t thirdHarmonicArray[2] = {harmonicsArray[4], harmonicsArray[5]}; // Array with (m,-m).
767  Int_t fourthHarmonicArray[2] = {harmonicsArray[6], harmonicsArray[7]}; // Array with (n,-n).
768 
769  Double_t tempoThirdHarmoNested = 0.;
770 
772  TComplex firstHarmonicComplex = (CalculateRecursion(2, firstHarmonicArray))/twoParticleDenominator; // Complex value of <2>_{k,-k}.
773  Double_t firstHarmonicValue = firstHarmonicComplex.Re(); // Value of <2>_{k,-k}.
774  fProfileCosineTwoParticles[0]->Fill(0.5, firstHarmonicValue, twoParticleDenominator);
775  fProfileCosineTwoParticles[0]->SetTitle(Form("#LT2#GT_{k,-k}, k = %d", firstHarmonicArray[0]));
776  fProfileCosineTwoParticles[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", firstHarmonicArray[0], firstHarmonicArray[1]));
777 
779  TComplex secondHarmonicComplex = (CalculateRecursion(2, secondHarmonicArray))/twoParticleDenominator; // Complex value of <2>_{l,-l}.
780  Double_t secondHarmonicValue = secondHarmonicComplex.Re(); // Value of <2>_{l,-l}.
781  fProfileCosineTwoParticles[1]->Fill(0.5, secondHarmonicValue, twoParticleDenominator);
782  fProfileCosineTwoParticles[1]->SetTitle(Form("#LT2#GT_{l,-l}, l = %d", secondHarmonicArray[0]));
783  fProfileCosineTwoParticles[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", secondHarmonicArray[0], secondHarmonicArray[1]));
784 
786  TComplex thirdHarmonicComplex = TComplex(0.,0.); // Complex value of <2>_{m,-m}.
787  Double_t thirdHarmonicValue = 0.; // Value of <2>_{m,-m}.
788  if (fNumberHarmonicsInSC >= 3)
789  {
790  thirdHarmonicComplex = (CalculateRecursion(2, thirdHarmonicArray))/twoParticleDenominator;
791  thirdHarmonicValue = thirdHarmonicComplex.Re();
792  fProfileCosineTwoParticles[2]->Fill(0.5, thirdHarmonicValue, twoParticleDenominator);
793  fProfileCosineTwoParticles[2]->SetTitle(Form("#LT2#GT_{m,-m}, m = %d", thirdHarmonicArray[0]));
794  fProfileCosineTwoParticles[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", thirdHarmonicArray[0], thirdHarmonicArray[1]));
795  } // End: if (fNumberHarmonicsInSC >= 3)
796 
798  TComplex fourthHarmonicComplex = TComplex(0.,0.); // Complex value of <2>_{n,-n}.
799  Double_t fourthHarmonicValue = 0.; // Value of <2>_{n,-n}.
800  if (fNumberHarmonicsInSC == 4)
801  {
802  fourthHarmonicComplex = (CalculateRecursion(2, fourthHarmonicArray))/twoParticleDenominator;
803  fourthHarmonicValue = fourthHarmonicComplex.Re();
804  fProfileCosineTwoParticles[3]->Fill(0.5, fourthHarmonicValue, twoParticleDenominator);
805  fProfileCosineTwoParticles[3]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", fourthHarmonicArray[0]));
806  fProfileCosineTwoParticles[3]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", fourthHarmonicArray[0], fourthHarmonicArray[1]));
807  } // End: if (fNumberHarmonicsInSC == 4)
808 
810  Int_t sumTwoHarmonicsArray[2] = {firstHarmonicArray[0]+secondHarmonicArray[0], firstHarmonicArray[1]+secondHarmonicArray[1]}; // Array with (k+l,-k-l).
811  TComplex sumTwoHarmonicsComplex = TComplex(0.,0.); // Complex value of <2>_{k+l,-k-l}.
812  Double_t sumTwoHarmonicsValue = 0.; // <2>_{k+l,-k-l}.
813  if (fNumberHarmonicsInSC == 2)
814  {
815  sumTwoHarmonicsComplex = (CalculateRecursion(2, sumTwoHarmonicsArray))/twoParticleDenominator;
816  sumTwoHarmonicsValue = sumTwoHarmonicsComplex.Re();
817  fProfileCosineTwoParticles[4]->Fill(0.5, sumTwoHarmonicsValue, twoParticleDenominator);
818  fProfileCosineTwoParticles[4]->SetTitle(Form("#LT2#GT_{k+l,-k-l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
819  fProfileCosineTwoParticles[4]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", sumTwoHarmonicsArray[0], sumTwoHarmonicsArray[1]));
820  } // End: if (fNumberHarmonicsInSC == 2)
821 
823  Int_t diffTwoHarmonicsArray[2] = {firstHarmonicArray[0]+secondHarmonicArray[1], firstHarmonicArray[1]+secondHarmonicArray[0]}; // Array with (k-l,-k+l).
824  TComplex diffTwoHarmonicsComplex = TComplex(0.,0.); // Complex value of <2>_{k-l,-k+l}.
825  Double_t diffTwoHarmonicsValue = 0.; // <2>_{k-l,-k+l}.
826  if (fNumberHarmonicsInSC == 2)
827  {
828  diffTwoHarmonicsComplex = (CalculateRecursion(2, diffTwoHarmonicsArray))/twoParticleDenominator;
829  diffTwoHarmonicsValue = diffTwoHarmonicsComplex.Re();
830  fProfileCosineTwoParticles[5]->Fill(0.5, diffTwoHarmonicsValue, twoParticleDenominator);
831  fProfileCosineTwoParticles[5]->SetTitle(Form("#LT2#GT_{k-l,-k+l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
832  fProfileCosineTwoParticles[5]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", diffTwoHarmonicsArray[0], diffTwoHarmonicsArray[1]));
833  } // End: if (fNumberHarmonicsInSC == 2)
834 
836  Double_t doubleCosineValueFirstCase = firstHarmonicValue*secondHarmonicValue;
837  fProfileTwoCosine[0]->Fill(0.5, doubleCosineValueFirstCase);
838  fProfileTwoCosine[0]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{l,-l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
839  fProfileTwoCosine[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], secondHarmonicArray[0], secondHarmonicArray[1]));
840 
842  Double_t doubleCosineValueSecondCase = 0.;
843  if (fNumberHarmonicsInSC >= 3)
844  {
845  doubleCosineValueSecondCase = firstHarmonicValue*thirdHarmonicValue;
846  fProfileTwoCosine[1]->Fill(0.5, doubleCosineValueSecondCase);
847  fProfileTwoCosine[1]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{m,-m}, k = %d, m = %d", firstHarmonicArray[0], thirdHarmonicArray[0]));
848  fProfileTwoCosine[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
849  } // End: if (fNumberHarmonicsInSC >= 3)
850 
852  Double_t doubleCosineValueThirdCase = 0.;
853  if (fNumberHarmonicsInSC >= 3)
854  {
855  doubleCosineValueThirdCase = secondHarmonicValue*thirdHarmonicValue;
856  fProfileTwoCosine[2]->Fill(0.5, doubleCosineValueThirdCase);
857  fProfileTwoCosine[2]->SetTitle(Form("#LT2#GT_{l,-l}#LT2#GT_{m,-m}, l = %d, m = %d", secondHarmonicArray[0], thirdHarmonicArray[0]));
858  fProfileTwoCosine[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", secondHarmonicArray[0], secondHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
859  } // End: if (fNumberHarmonicsInSC >= 3)
860 
862  Double_t doubleCosineValueFourthCase = 0.;
863  if (fNumberHarmonicsInSC >= 3)
864  {
865  doubleCosineValueFourthCase = firstHarmonicValue*secondHarmonicValue*thirdHarmonicValue;
866  fProfileTwoCosine[3]->Fill(0.5, doubleCosineValueFourthCase);
867  fProfileTwoCosine[3]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{l,-l}#LT2#GT_{m,-m}, k = %d, l = %d, m = %d", firstHarmonicArray[0], secondHarmonicArray[0], thirdHarmonicArray[0]));
868  fProfileTwoCosine[3]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], secondHarmonicArray[0], secondHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
869  } // End: if (fNumberHarmonicsInSC >= 3)
870 
873  {
874  // 1st harmonic: <2>_{k,-k}.
875  Double_t firstHarmonicValueNested = ComputeTwoNestedLoops(nParticles, firstHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[0]);
876  fProfileCosineTwoNestedLoops[0]->SetTitle(Form("#LT2#GT_{k,-k}, k = %d", firstHarmonicArray[0]));
877  fProfileCosineTwoNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", firstHarmonicArray[0], firstHarmonicArray[1]));
878 
879  // 2nd harmonic: <2>_{l,-l}.
880  Double_t secondHarmonicValueNested = ComputeTwoNestedLoops(nParticles, secondHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[1]);
881  fProfileCosineTwoNestedLoops[1]->SetTitle(Form("#LT2#GT_{l,-l}, l = %d", secondHarmonicArray[0]));
882  fProfileCosineTwoNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", secondHarmonicArray[0], secondHarmonicArray[1]));
883 
884  // 3rd harmonic: <2>_{m,-m}.
885  Double_t thirdHarmonicValueNested = 0.;
886  if (fNumberHarmonicsInSC >= 3)
887  {
888  thirdHarmonicValueNested = ComputeTwoNestedLoops(nParticles, thirdHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[2]);
889  fProfileCosineTwoNestedLoops[2]->SetTitle(Form("#LT2#GT_{m,-m}, m = %d", thirdHarmonicArray[0]));
890  fProfileCosineTwoNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", thirdHarmonicArray[0], thirdHarmonicArray[1]));
891  } // End: if (fNumberHarmonicsInSC >= 3)
892  tempoThirdHarmoNested = thirdHarmonicValueNested;
893 
894  // 4th harmonic: <2>_{n,-n}.
895  Double_t fourthHarmonicValueNested = 0.;
896  fourthHarmonicValueNested++; // Dummy line since this variable is not useful yet.
897  if (fNumberHarmonicsInSC == 4)
898  {
899  fourthHarmonicValueNested = ComputeTwoNestedLoops(nParticles, fourthHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[3]);
900  fProfileCosineTwoNestedLoops[3]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", fourthHarmonicArray[0]));
901  fProfileCosineTwoNestedLoops[3]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", fourthHarmonicArray[0], fourthHarmonicArray[1]));
902  } // End: if (fNumberHarmonicsInSC == 4)
903 
904  // Sum first two harmonics: <2>_{k+l,-k-l}.
905  Double_t sumTwoHarmonicsValueNested = 0.;
906  sumTwoHarmonicsValueNested++; // Dummy line since this variable is not useful yet.
907  if (fNumberHarmonicsInSC == 2)
908  {
909  sumTwoHarmonicsValueNested = ComputeTwoNestedLoops(nParticles, sumTwoHarmonicsArray, angles, weights, fProfileCosineTwoNestedLoops[4]);
910  fProfileCosineTwoNestedLoops[4]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", sumTwoHarmonicsArray[0]));
911  fProfileCosineTwoNestedLoops[4]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", sumTwoHarmonicsArray[0], sumTwoHarmonicsArray[1]));
912  } // End: if (fNumberHarmonicsInSC == 2)
913 
914  // Difference first two harmonics: <2>_{k-l,-k+l}.
915  Double_t diffTwoHarmonicsValueNested = 0.;
916  diffTwoHarmonicsValueNested++; // Dummy line since this variable is not useful yet.
917  if (fNumberHarmonicsInSC == 2)
918  {
919  diffTwoHarmonicsValueNested = ComputeTwoNestedLoops(nParticles, diffTwoHarmonicsArray, angles, weights, fProfileCosineTwoNestedLoops[5]);
920  fProfileCosineTwoNestedLoops[5]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", diffTwoHarmonicsArray[0]));
921  fProfileCosineTwoNestedLoops[5]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", diffTwoHarmonicsArray[0], diffTwoHarmonicsArray[1]));
922  } // End: if (fNumberHarmonicsInSC == 2)
923 
924  // <<2>_{k,-k}<2>_{l,-l}>.
925  Double_t doubleCosineFirstCaseNested = firstHarmonicValueNested*secondHarmonicValueNested;
926  fProfileTwoCosineNestedLoops[0]->Fill(0.5, doubleCosineFirstCaseNested);
927  fProfileTwoCosineNestedLoops[0]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{l,-l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
928  fProfileTwoCosineNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], secondHarmonicArray[0], secondHarmonicArray[1]));
929 
930  // <2>_{k,-k}<2>_{m,-m}.
931  Double_t doubleCosineSecondCaseNested = 0.;
932  if (fNumberHarmonicsInSC >= 3)
933  {
934  doubleCosineSecondCaseNested = firstHarmonicValueNested*thirdHarmonicValueNested;
935  fProfileTwoCosineNestedLoops[1]->Fill(0.5, doubleCosineSecondCaseNested);
936  fProfileTwoCosineNestedLoops[1]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{m,-m}, k = %d, m = %d", firstHarmonicArray[0], thirdHarmonicArray[0]));
937  fProfileTwoCosineNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
938  } // End: if (fNumberHarmonicsInSC >= 3)
939 
940  // <2>_{l,-l}<2>_{m,-m}.
941  Double_t doubleCosineThirdCaseNested = 0.;
942  if (fNumberHarmonicsInSC >= 3)
943  {
944  doubleCosineThirdCaseNested = secondHarmonicValueNested*thirdHarmonicValueNested;
945  fProfileTwoCosineNestedLoops[2]->Fill(0.5, doubleCosineThirdCaseNested);
946  fProfileTwoCosineNestedLoops[2]->SetTitle(Form("#LT2#GT_{l,-l}#LT2#GT_{m,-m}, l = %d, m = %d", secondHarmonicArray[0], thirdHarmonicArray[0]));
947  fProfileTwoCosineNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", secondHarmonicArray[0], secondHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
948  } // End: if (fNumberHarmonicsInSC >= 3)
949 
950  // <2>_{k,-k}<2>_{l,-l}<2>_{m,-m}.
951  Double_t doubleCosineFourthCaseNested = 0.;
952  if (fNumberHarmonicsInSC >= 3)
953  {
954  doubleCosineFourthCaseNested = firstHarmonicValueNested*secondHarmonicValueNested*thirdHarmonicValueNested;
955  fProfileTwoCosineNestedLoops[3]->Fill(0.5, doubleCosineFourthCaseNested);
956  fProfileTwoCosineNestedLoops[3]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{l,-l}#LT2#GT_{m,-m}, k = %d, l = %d, m = %d", firstHarmonicArray[0], secondHarmonicArray[0], thirdHarmonicArray[0]));
957  fProfileTwoCosineNestedLoops[3]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], secondHarmonicArray[0], secondHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
958  } // End: if (fNumberHarmonicsInSC >= 3)
959 
960  // Reset of the observables for the two nested loops.
961  firstHarmonicValueNested = 0.;
962  secondHarmonicValueNested = 0.;
963  thirdHarmonicValueNested = 0.;
964  fourthHarmonicValueNested = 0.;
965  sumTwoHarmonicsValueNested = 0.;
966  diffTwoHarmonicsValueNested = 0.;
967  doubleCosineFirstCaseNested = 0.;
968  doubleCosineSecondCaseNested = 0.;
969  doubleCosineThirdCaseNested = 0.;
970  doubleCosineFourthCaseNested = 0.;
971  } // End: if (fComputeNestedLoops)
972 
973 //**************************************************************************************//
974 // Compute the 3-particle correlations.
976  Int_t threeParticleCaseOneArray[3] = {harmonicsArray[0]+harmonicsArray[2], harmonicsArray[1], harmonicsArray[3]};
977  TComplex threeParticleCaseOneComplex = TComplex(0.,0.); // Complex value of <3>_{k+l,-k,-l}.
978  Double_t threeParticleCaseOneValue = 0.; // Value of <3>_{k+l,-k,-l}.
979  if (fNumberHarmonicsInSC == 2)
980  {
981  threeParticleCaseOneComplex = (CalculateRecursion(3, threeParticleCaseOneArray))/threeParticleDenominator;
982  threeParticleCaseOneValue = threeParticleCaseOneComplex.Re();
983  fProfileCosineThreeParticles[0]->Fill(0.5, threeParticleCaseOneValue, threeParticleDenominator);
984  fProfileCosineThreeParticles[0]->SetTitle(Form("#LT3#GT_{k+l,-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
985  fProfileCosineThreeParticles[0]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseOneArray[0], threeParticleCaseOneArray[1], threeParticleCaseOneArray[2]));
986  } // End: if (fNumberHarmonicsInSC == 2)
987 
989  Int_t threeParticleCaseTwoArray[3] = {harmonicsArray[0]+harmonicsArray[3], harmonicsArray[1], harmonicsArray[2]};
990  TComplex threeParticleCaseTwoComplex = TComplex(0.,0.); // Complex value of <3>_{k+l,-k,-l}.
991  Double_t threeParticleCaseTwoValue = 0.; // Value of <3>_{k-l,-k,+l}.
992  if (fNumberHarmonicsInSC == 2)
993  {
994  threeParticleCaseTwoComplex = (CalculateRecursion(3, threeParticleCaseTwoArray))/threeParticleDenominator;
995  threeParticleCaseTwoValue = threeParticleCaseTwoComplex.Re();
996  fProfileCosineThreeParticles[1]->Fill(0.5, threeParticleCaseTwoValue, threeParticleDenominator);
997  fProfileCosineThreeParticles[1]->SetTitle(Form("#LT3#GT_{k-l,-k,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
998  fProfileCosineThreeParticles[1]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseTwoArray[0], threeParticleCaseTwoArray[1], threeParticleCaseTwoArray[2]));
999  } // End: if (fNumberHarmonicsInSC == 2)
1000 
1002  Int_t threeParticleCaseThreeArray[3] = {harmonicsArray[0], harmonicsArray[2]+harmonicsArray[1], harmonicsArray[3]};
1003  TComplex threeParticleCaseThreeComplex = TComplex(0.,0.); // Complex value of <3>_{k,l-k,-l}.
1004  Double_t threeParticleCaseThreeValue = 0.; // Value of <3>_{k,l-k,-l}.
1005  if (fNumberHarmonicsInSC == 2)
1006  {
1007  threeParticleCaseThreeComplex = (CalculateRecursion(3, threeParticleCaseThreeArray))/threeParticleDenominator;
1008  threeParticleCaseThreeValue = threeParticleCaseThreeComplex.Re();
1009  fProfileCosineThreeParticles[2]->Fill(0.5, threeParticleCaseThreeValue, threeParticleDenominator);
1010  fProfileCosineThreeParticles[2]->SetTitle(Form("#LT3#GT_{k,l-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1011  fProfileCosineThreeParticles[2]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseThreeArray[0], threeParticleCaseThreeArray[1], threeParticleCaseThreeArray[2]));
1012  } // End: if (fNumberHarmonicsInSC == 2)
1013 
1015  Int_t threeParticleCaseFourArray[3] = {harmonicsArray[0], harmonicsArray[1]+harmonicsArray[3], harmonicsArray[2]};
1016  TComplex threeParticleCaseFourComplex = TComplex(0.,0.); // Complex value of <3>_{k,-k-l,l}.
1017  Double_t threeParticleCaseFourValue = 0.; // Value of <3>_{k,-k-l,l}.
1018  if (fNumberHarmonicsInSC == 2)
1019  {
1020  threeParticleCaseFourComplex = (CalculateRecursion(3, threeParticleCaseFourArray))/threeParticleDenominator;
1021  threeParticleCaseFourValue = threeParticleCaseFourComplex.Re();
1022  fProfileCosineThreeParticles[3]->Fill(0.5, threeParticleCaseFourValue, threeParticleDenominator);
1023  fProfileCosineThreeParticles[3]->SetTitle(Form("#LT3#GT_{k,-k-l,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1024  fProfileCosineThreeParticles[3]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseFourArray[0], threeParticleCaseFourArray[1], threeParticleCaseFourArray[2]));
1025  } // End: if (fNumberHarmonicsInSC == 2)
1026 
1028  if (fComputeNestedLoops)
1029  {
1030  // Case one: <3>_{k+l,-k,-l}.
1031  Double_t threeParticleCaseOneNested = 0.;
1032  threeParticleCaseOneNested++; // Dummy line since this variable is not useful yet.
1033  if (fNumberHarmonicsInSC == 2)
1034  {
1035  threeParticleCaseOneNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseOneArray, angles, weights, fProfileCosineThreeNestedLoops[0]);
1036  fProfileCosineThreeNestedLoops[0]->SetTitle(Form("#LT3#GT_{k+l,-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1037  fProfileCosineThreeNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseOneArray[0], threeParticleCaseOneArray[1], threeParticleCaseOneArray[2]));
1038  } // End: if (fNumberHarmonicsInSC == 2)
1039 
1040  // Case two: <3>_{k-l,-k,+l}.
1041  Double_t threeParticleCaseTwoNested = 0.;
1042  threeParticleCaseTwoNested++; // Dummy line since this variable is not useful yet.
1043  if (fNumberHarmonicsInSC == 2)
1044  {
1045  threeParticleCaseTwoNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseTwoArray, angles, weights, fProfileCosineThreeNestedLoops[1]);
1046  fProfileCosineThreeNestedLoops[1]->SetTitle(Form("#LT3#GT_{k-l,-k,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1047  fProfileCosineThreeNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseTwoArray[0], threeParticleCaseTwoArray[1], threeParticleCaseTwoArray[2]));
1048  } // End: if (fNumberHarmonicsInSC == 2)
1049 
1050  // Case three: <3>_{k,l-k,-l}.
1051  Double_t threeParticleCaseThreeNested = 0.;
1052  threeParticleCaseThreeNested++; // Dummy line since this variable is not useful yet.
1053  if (fNumberHarmonicsInSC == 2)
1054  {
1055  threeParticleCaseThreeNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseThreeArray, angles, weights, fProfileCosineThreeNestedLoops[2]);
1056  fProfileCosineThreeNestedLoops[2]->SetTitle(Form("#LT3#GT_{k,l-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1057  fProfileCosineThreeNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseThreeArray[0], threeParticleCaseThreeArray[1], threeParticleCaseThreeArray[2]));
1058  } // End: if (fNumberHarmonicsInSC == 2)
1059 
1060  // Case four: <3>_{k,-k-l,l}.
1061  Double_t threeParticleCaseFourNested = 0.;
1062  threeParticleCaseFourNested++; // Dummy line since this variable is not useful yet.
1063  if (fNumberHarmonicsInSC == 2)
1064  {
1065  threeParticleCaseFourNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseFourArray, angles, weights, fProfileCosineThreeNestedLoops[3]);
1066  fProfileCosineThreeNestedLoops[3]->SetTitle(Form("#LT3#GT_{k,-k-l,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1067  fProfileCosineThreeNestedLoops[3]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseFourArray[0], threeParticleCaseFourArray[1], threeParticleCaseFourArray[2]));
1068  } // End: if (fNumberHarmonicsInSC == 2)
1069 
1070  // Reset of the observables for the three nested loops.
1071  threeParticleCaseOneNested = 0.;
1072  threeParticleCaseTwoNested = 0.;
1073  threeParticleCaseThreeNested = 0.;
1074  threeParticleCaseFourNested = 0.;
1075  } // End: if (fComputeNestedLoops)
1076 
1077 //**************************************************************************************//
1078 // Compute the 4-particle correlations.
1080  Int_t fourParticleCaseOneArray[4] = {harmonicsArray[0], harmonicsArray[1], harmonicsArray[2], harmonicsArray[3]};
1081  TComplex fourParticleCaseOneComplex = (CalculateRecursion(4, fourParticleCaseOneArray))/fourParticleDenominator; // Complex value of <4>_{k,-k,l,-l}.
1082  Double_t fourParticleCaseOneValue = fourParticleCaseOneComplex.Re(); // Value of <4>_{k,-k,l,-l}.
1083  fProfileCosineFourParticles[0]->Fill(0.5, fourParticleCaseOneValue, fourParticleDenominator);
1084  fProfileCosineFourParticles[0]->SetTitle(Form("#LT4#GT_{k,-k,l,-l}, k = %d, l = %d", fourParticleCaseOneArray[0], fourParticleCaseOneArray[2]));
1085  fProfileCosineFourParticles[0]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseOneArray[0], fourParticleCaseOneArray[1],fourParticleCaseOneArray[2],fourParticleCaseOneArray[3]));
1086 
1088  Int_t fourParticleCaseTwoArray[4] = {harmonicsArray[0], harmonicsArray[1], harmonicsArray[4], harmonicsArray[5]};
1089  TComplex fourParticleCaseTwoComplex = TComplex(0.,0.); // Complex value of <4>_{k,-k,m,-m}.
1090  Double_t fourParticleCaseTwoValue = 0.; // Value of <4>_{k,-k,m,-m}.
1091  if (fNumberHarmonicsInSC >= 3)
1092  {
1093  fourParticleCaseTwoComplex = (CalculateRecursion(4, fourParticleCaseTwoArray))/fourParticleDenominator;
1094  fourParticleCaseTwoValue = fourParticleCaseTwoComplex.Re();
1095  fProfileCosineFourParticles[1]->Fill(0.5, fourParticleCaseTwoValue, fourParticleDenominator);
1096  fProfileCosineFourParticles[1]->SetTitle(Form("#LT4#GT_{k,-k,m,-m}, k = %d, m = %d", fourParticleCaseTwoArray[0], fourParticleCaseTwoArray[2]));
1097  fProfileCosineFourParticles[1]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseTwoArray[0], fourParticleCaseTwoArray[1],fourParticleCaseTwoArray[2],fourParticleCaseTwoArray[3]));
1098  } // End: if (fNumberHarmonicsInSC >= 3)
1099 
1101  Int_t fourParticleCaseThreeArray[4] = {harmonicsArray[2], harmonicsArray[3], harmonicsArray[4], harmonicsArray[5]};
1102  TComplex fourParticleCaseThreeComplex = TComplex(0.,0.); // Complex value of <4>_{l,-l,m,-m}.
1103  Double_t fourParticleCaseThreeValue = 0.; // Value of <4>_{l,-l,m,-m}.
1104  if (fNumberHarmonicsInSC >= 3)
1105  {
1106  fourParticleCaseThreeComplex = (CalculateRecursion(4, fourParticleCaseThreeArray))/fourParticleDenominator;
1107  fourParticleCaseThreeValue = fourParticleCaseThreeComplex.Re();
1108  fProfileCosineFourParticles[2]->Fill(0.5, fourParticleCaseThreeValue, fourParticleDenominator);
1109  fProfileCosineFourParticles[2]->SetTitle(Form("#LT4#GT_{l,-l,m,-m}, l = %d, m = %d", fourParticleCaseThreeArray[0], fourParticleCaseThreeArray[2]));
1110  fProfileCosineFourParticles[2]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseThreeArray[0], fourParticleCaseThreeArray[1],fourParticleCaseThreeArray[2],fourParticleCaseThreeArray[3]));
1111  } // End: if (fNumberHarmonicsInSC >= 3)
1112 
1114 
1116  Double_t doubleCosineValueFifthCase = 0.;
1117  if (fNumberHarmonicsInSC >= 3)
1118  {
1119  doubleCosineValueFifthCase = fourParticleCaseOneValue*thirdHarmonicValue;
1120  fProfileTwoCosine[4]->Fill(0.5, doubleCosineValueFifthCase);
1121  fProfileTwoCosine[4]->SetTitle(Form("#LT4#GT_{k,-k,l,-l}#LT2#GT_{m,-m}, k = %d, l = %d, m = %d", fourParticleCaseOneArray[0], fourParticleCaseOneArray[2], thirdHarmonicArray[0]));
1122  fProfileTwoCosine[4]->GetYaxis()->SetTitle(Form("#LT#LT4#GT_{%d,%d,%d,%d}#LT2#GT_{%d,%d}#GT", fourParticleCaseOneArray[0], fourParticleCaseOneArray[1], fourParticleCaseOneArray[2], fourParticleCaseOneArray[3], thirdHarmonicArray[0], thirdHarmonicArray[1]));
1123  } // End: if (fNumberHarmonicsInSC >= 3)
1124 
1126  if (fComputeNestedLoops)
1127  {
1128  // Case one: <4>_{k,-k,l,-l}.
1129  Double_t fourParticleCaseOneNested = ComputeFourNestedLoops(nParticles, fourParticleCaseOneArray, angles, weights, fProfileCosineFourNestedLoops[0]);
1130  fProfileCosineFourNestedLoops[0]->SetTitle(Form("#LT4#GT_{k,-k,l,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1131  fProfileCosineFourNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseOneArray[0], fourParticleCaseOneArray[1], fourParticleCaseOneArray[2],fourParticleCaseOneArray[3]));
1132 
1133  // Case two: <4>_{k,m,-k,-m}.
1134  Double_t fourParticleCaseTwoNested = 0.;
1135  fourParticleCaseTwoNested++; // Dummy line since this variable is not useful yet.
1136  if (fNumberHarmonicsInSC >= 3)
1137  {
1138  fourParticleCaseTwoNested = ComputeFourNestedLoops(nParticles, fourParticleCaseTwoArray, angles, weights, fProfileCosineFourNestedLoops[1]);
1139  fProfileCosineFourNestedLoops[1]->SetTitle(Form("#LT4#GT_{k,-k,m,-m}, k = %d, m = %d", harmonicsArray[0], harmonicsArray[4]));
1140  fProfileCosineFourNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseTwoArray[0], fourParticleCaseTwoArray[1], fourParticleCaseTwoArray[2],fourParticleCaseTwoArray[3]));
1141  } // End: if (fNumberHarmonicsInSC >= 3)
1142 
1143  // Case three: <4>_{l,m,-l,-m}.
1144  Double_t fourParticleCaseThreeNested = 0.;
1145  fourParticleCaseThreeNested++; // Dummy line since this variable is not useful yet.
1146  if (fNumberHarmonicsInSC >= 3)
1147  {
1148  fourParticleCaseThreeNested = ComputeFourNestedLoops(nParticles, fourParticleCaseThreeArray, angles, weights, fProfileCosineFourNestedLoops[2]);
1149  fProfileCosineFourNestedLoops[2]->SetTitle(Form("#LT4#GT_{l,-l,m,-m}, l = %d, m = %d", harmonicsArray[2], harmonicsArray[4]));
1150  fProfileCosineFourNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseThreeArray[0], fourParticleCaseThreeArray[1], fourParticleCaseThreeArray[2],fourParticleCaseThreeArray[3]));
1151  } // End: if (fNumberHarmonicsInSC >= 3)
1152 
1153  // TBI: cases with four harmonics...
1154 
1155  // <4>_{k,-k,l,-l}<2>_{m,-m}.
1156  Double_t doubleCosineFifthCaseNested = 0.;
1157  doubleCosineFifthCaseNested++; // Dummy line since this variable is not useful yet.
1158  if (fNumberHarmonicsInSC >= 3)
1159  {
1160  doubleCosineFifthCaseNested = fourParticleCaseOneNested*tempoThirdHarmoNested;
1161  fProfileTwoCosineNestedLoops[4]->Fill(0.5, doubleCosineFifthCaseNested);
1162  fProfileTwoCosineNestedLoops[4]->SetTitle(Form("#LT4#GT_{k,-k,l,-l}#LT2#GT_{m,-m}, k = %d, l = %d, m = %d", fourParticleCaseOneArray[0], fourParticleCaseOneArray[2], thirdHarmonicArray[0]));
1163  fProfileTwoCosineNestedLoops[4]->GetYaxis()->SetTitle(Form("#LT#LT4#GT_{%d,%d,%d,%d}#LT2#GT_{%d,%d}#GT", fourParticleCaseOneArray[0], fourParticleCaseOneArray[1], fourParticleCaseOneArray[2], fourParticleCaseOneArray[3], thirdHarmonicArray[0], thirdHarmonicArray[1]));
1164  } // End: if (fNumberHarmonicsInSC >= 3)
1165 
1166  // Reset of the observables for the four nested loops.
1167  fourParticleCaseOneNested = 0.;
1168  fourParticleCaseTwoNested = 0.;
1169  fourParticleCaseThreeNested = 0.;
1170  //
1171  doubleCosineFifthCaseNested = 0.;
1172  } // End: if (fComputeNestedLoops)
1173 
1174 //**************************************************************************************//
1175 // Compute the 6-particle correlations.
1177  Int_t sixParticleCaseOneArray[6] = {harmonicsArray[0], harmonicsArray[1], harmonicsArray[2], harmonicsArray[3], harmonicsArray[4], harmonicsArray[5]};
1178  TComplex sixParticleCaseOneComplex = TComplex(0.,0.); // Complex value of <6>_{k,-k,l,-l,m,-m}.
1179  Double_t sixParticleCaseOneValue = 0.; // Value of <6>_{k,-k,l,-l,m,-m}.
1180  if (fNumberHarmonicsInSC >= 3)
1181  {
1182  sixParticleCaseOneComplex = (CalculateRecursion(6, sixParticleCaseOneArray))/sixParticleDenominator;
1183  sixParticleCaseOneValue = sixParticleCaseOneComplex.Re();
1184  fProfileCosineSixParticles[0]->Fill(0.5, sixParticleCaseOneValue, sixParticleDenominator);
1185  fProfileCosineSixParticles[0]->SetTitle(Form("#LT6#GT_{k,-k,l,-l,m,-m}, k = %d, l = %d, m = %d", sixParticleCaseOneArray[0], sixParticleCaseOneArray[2], sixParticleCaseOneArray[4]));
1186  fProfileCosineSixParticles[0]->GetYaxis()->SetTitle(Form("#LT#LT6#GT#GT_{%d,%d,%d,%d,%d,%d}", sixParticleCaseOneArray[0], sixParticleCaseOneArray[1], sixParticleCaseOneArray[2], sixParticleCaseOneArray[3], sixParticleCaseOneArray[4], sixParticleCaseOneArray[5]));
1187  } // End: if (fNumberHarmonicsInSC >= 3)
1188 
1189 //**************************************************************************************//
1190 // Compute the 8-particle correlations.
1192 
1193 //**************************************************************************************//
1194 // Reset of the variables to zero.
1195  twoParticleDenominator = 0.;
1196  threeParticleDenominator = 0.;
1197  fourParticleDenominator = 0.;
1198  sixParticleDenominator = 0.;
1199  eightParticleDenominator = 0.;
1200  tempoThirdHarmoNested = 0.;
1201 
1202  firstHarmonicComplex = TComplex(0.,0.);
1203  firstHarmonicValue = 0.;
1204  secondHarmonicComplex = TComplex(0.,0.);
1205  secondHarmonicValue = 0.;
1206  thirdHarmonicComplex = TComplex(0.,0.);
1207  thirdHarmonicValue = 0.;
1208  fourthHarmonicComplex = TComplex(0.,0.);
1209  fourthHarmonicValue = 0.;
1210 
1211  doubleCosineValueFirstCase = 0.,
1212  doubleCosineValueSecondCase = 0.;
1213  doubleCosineValueThirdCase = 0.;
1214 
1215  threeParticleCaseOneComplex = TComplex(0.,0.);
1216  threeParticleCaseOneValue = 0.;
1217  threeParticleCaseTwoComplex = TComplex(0.,0.);
1218  threeParticleCaseTwoValue = 0.;
1219  threeParticleCaseThreeComplex = TComplex(0.,0.);
1220  threeParticleCaseThreeValue = 0.;
1221  threeParticleCaseFourComplex = TComplex(0.,0.);
1222  threeParticleCaseFourValue = 0.;
1223 
1224  fourParticleCaseOneComplex = TComplex(0.,0.);
1225  fourParticleCaseOneValue = 0.;
1226  fourParticleCaseTwoComplex = TComplex(0.,0.);
1227  fourParticleCaseTwoValue = 0.;
1228  fourParticleCaseThreeComplex = TComplex(0.,0.);
1229  fourParticleCaseThreeValue = 0.;
1230  // TBI: cases with four possible harmonics: kn, ln, mn...
1231  doubleCosineValueFifthCase = 0.;
1232 
1233  sixParticleCaseOneComplex = TComplex(0.,0.);
1234  sixParticleCaseOneValue = 0.;
1235 } // End: GSCfullAnalysis(long long, Double_t[], Double_t[])
1236 
1237 //--------------------------------------------------------------------------------------//
1238 
1240 {
1241 // Simplify the use of the Q-vectors.
1242  if (n >= 0) {return fQvectors[n][p];}
1243  return TComplex::Conjugate(fQvectors[-n][p]);
1244 } // End: Q(Int_t, Int_t)
1245 
1246 //--------------------------------------------------------------------------------------//
1247 
1249 {
1250 // Calculate multi-particle correlators by using recursion (an improved faster version) originally developed by Kristjan Gulbrandsen (gulbrand@nbi.dk).
1251  Int_t nm1 = n-1;
1252  TComplex c(Q(harmonic[nm1], mult));
1253  if (nm1 == 0) return c;
1254  c *= CalculateRecursion(nm1, harmonic);
1255  if (nm1 == skip) return c;
1256 
1257  Int_t multp1 = mult+1;
1258  Int_t nm2 = n-2;
1259  Int_t counter1 = 0;
1260  Int_t hhold = harmonic[counter1];
1261  harmonic[counter1] = harmonic[nm2];
1262  harmonic[nm2] = hhold + harmonic[nm1];
1263  TComplex c2(CalculateRecursion(nm1, harmonic, multp1, nm2));
1264  Int_t counter2 = n-3;
1265 
1266  while (counter2 >= skip) {
1267  harmonic[nm2] = harmonic[counter1];
1268  harmonic[counter1] = hhold;
1269  ++counter1;
1270  hhold = harmonic[counter1];
1271  harmonic[counter1] = harmonic[nm2];
1272  harmonic[nm2] = hhold + harmonic[nm1];
1273  c2 += CalculateRecursion(nm1, harmonic, multp1, counter2);
1274  --counter2;
1275  }
1276  harmonic[nm2] = harmonic[counter1];
1277  harmonic[counter1] = hhold;
1278 
1279  if (mult == 1) return c-c2;
1280  return c-Double_t(mult)*c2;
1281 } // End: CalculateRecursion(Int_t, Int_t[], Int_t, Int_t)
1282 
1283 //--------------------------------------------------------------------------------------//
1284 
1285 Double_t AliAnalysisTaskTwoMultiCorrelations::ComputeTwoNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
1286 {
1287 // Calculate the two-particle correlations using two nested loops.
1288  Double_t twoParticleCosine = 0.; // cos(k(phi_1 - phi_2)).
1289  Double_t twoParticleWeight = 1.; // Total particle weights.
1290  TProfile twoProfile; // Returns <cos(k(phi1-phi2))>.
1291  twoProfile.Sumw2();
1292 
1293  for (long long m = 0; m < nParticles; m++)
1294  {
1295  for (long long n = 0; n < nParticles; n++)
1296  {
1297  // Remove the autocorrelations m==n.
1298  if (m == n) {continue;}
1299  // Compute cos(k(phi_1-phi_2)).
1300  twoParticleCosine = TMath::Cos(harmonic[0]*angles[m] + harmonic[1]*angles[n]);
1301  twoParticleWeight = weights[m] + weights[n];
1302  profile->Fill(0.5, twoParticleCosine, twoParticleWeight);
1303  twoProfile.Fill(0.5, twoParticleCosine, twoParticleWeight);
1304 
1305  // Reset the variables to default.
1306  twoParticleCosine = 0.;
1307  twoParticleWeight = 1.;
1308  } // End: for (long long n = 0; n < nParticles; n++)
1309  } // End: for (long long m = 0; m < nParticles; m++)
1310 
1311 // Return <cos(k(phi_1-phi_2))>.
1312  return twoProfile.GetBinContent(1);
1313 } // End: ComputeTwoNestedLoops(long long, Int_t*, Double_t[], Double_t[], TProfile*)
1314 
1315 //--------------------------------------------------------------------------------------//
1316 
1317 Double_t AliAnalysisTaskTwoMultiCorrelations::ComputeThreeNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
1318 {
1319 // Calculate the three-particle correlations using three nested loops.
1320  Double_t threeParticleCosine = 0.; // cos(kphi_1 +lphi_2 +mphi_3).
1321  Double_t threeParticleWeight = 1.; // Total particle weights.
1322  TProfile threeProfile; // Returns <cos(kphi_1 +lphi_2 +mphi_3)>.
1323  threeProfile.Sumw2();
1324 
1325  for (long long k = 0; k < nParticles; k++)
1326  {
1327  for (long long l = 0; l < nParticles; l++)
1328  {
1329  if (k == l) {continue;}
1330 
1331  for (long long m = 0; m < nParticles; m++)
1332  {
1333  if ((k == m) || (l == m)) {continue;}
1334  // Compute cos(kphi_1 +lphi_2 +mphi_3).
1335  threeParticleCosine = TMath::Cos(harmonic[0]*angles[k] + harmonic[1]*angles[l] + harmonic[2]*angles[m]);
1336  threeParticleWeight = weights[k] + weights[l] + weights[m];
1337  profile->Fill(0.5, threeParticleCosine, threeParticleWeight);
1338  threeProfile.Fill(0.5, threeParticleCosine, threeParticleWeight);
1339 
1340  // Reset the variables to default.
1341  threeParticleCosine = 0.;
1342  threeParticleWeight = 1.;
1343  } // End: for (long long m = 0; m < nParticles; m++)
1344  } // End: for (long long l = 0; l < nParticles; l++)
1345 } // End: for (long long k = 0; k < nParticles; k++)
1346 
1347  // Return <cos(kphi_1 +lphi_2 +mphi_3)>.
1348  return threeProfile.GetBinContent(1);
1349 } // End: ComputeThreeNestedLoops(long long, Int_t*, Double_t[], Double_t[], TProfile*)
1350 
1351 //--------------------------------------------------------------------------------------//
1352 
1353 Double_t AliAnalysisTaskTwoMultiCorrelations::ComputeFourNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
1354 {
1355 // Calculate the four-particle correlations using four nested loops.
1356  Double_t fourParticleCosine = 0.; // cos(kphi_1 +lphi_2 +mphi_3 +nphi_4).
1357  Double_t fourParticleWeight = 1.; // Total particle weights.
1358  TProfile fourProfile; // Returns <cos(kphi_1 +lphi_2 +mphi_3 +nphi_4)>.
1359  fourProfile.Sumw2();
1360 
1361  for (long long k = 0; k < nParticles; k++)
1362  {
1363  for (long long l = 0; l < nParticles; l++)
1364  {
1365  if (k == l) {continue;}
1366  for (long long m = 0; m < nParticles; m++)
1367  {
1368  if ((k == m) || (l == m)) {continue;}
1369  for (long long n = 0; n < nParticles; n++)
1370  {
1371  if ((k == n) || (l == n) || (m == n)) {continue;}
1372  // Computate cos(k(phi_1 - phi_2) + l(phi_3-phi_4)).
1373  fourParticleCosine = TMath::Cos(harmonic[0]*angles[k] + harmonic[1]*angles[l] + harmonic[2]*angles[m] + harmonic[3]*angles[n]);
1374  fourParticleWeight = weights[k] + weights[l] + weights[m] + weights[n];
1375  profile->Fill(0.5, fourParticleCosine, fourParticleWeight);
1376  fourProfile.Fill(0.5, fourParticleCosine, fourParticleWeight);
1377 
1378  // Reset the variables to default.
1379  fourParticleCosine = 0.;
1380  fourParticleWeight = 1.;
1381  } // End: for (long long n = 0; n < nParticles; n++)
1382  } // End: for (long long m = 0; m < nParticles; m++)
1383  } // End: for (long long l = 0; l < nParticles; l++)
1384  } // End: for (long long k = 0; k < nParticles; k++)
1385 
1386  // Return <cos(kphi_1 +lphi_2 +mphi_3 +nphi_4)>.
1387  return fourProfile.GetBinContent(1);
1388 } // End: Double_t ComputeFourNestedLoops
1389 
1390 //======================================================================================//
1391 // Methods called in Terminate(Option_t *).
1392 //--------------------------------------------------------------------------------------//
1393 
1394 //======================================================================================//
1395 
TProfile * fProfileCosineTwoNestedLoops[6]
<2>_{j,-j}, j: k,l,m,n,(k+l),(k-l).
TProfile * fProfileCosineTwoParticles[6]
Azimuthal angles distribution.
double Double_t
Definition: External.C:58
TProfile * fProfileCosineThreeParticles[4]
N<<2>_{i,-i}<2>_{j,-j}> with nested loops.
TProfile * fProfileCosineSixParticles[4]
<4>_{i,j,-i,-j} with nested loops.
TCanvas * c
Definition: TestFitELoss.C:172
virtual void CalculateQvectors(long long nParticles, Double_t angles[], Double_t weights[])
TH1D * fHistoPhiPostCuts
Pseudorapidity distribution.
int Int_t
Definition: External.C:63
TProfile * fProfileCosineEightParticles
<6>_{h,i,j,-h,-i,-j}, hij: klm,kln,kmn,lmn.
Definition: External.C:212
TH1D * fHistoEtaPreCuts
Transverse momentum distribution.
TH1D * fHistoMultiplicityPreCuts
Azimuthal angles distribution.
TH1D * fHistoMultiplicityPostCuts
Multiplicity distribution.
TProfile * fProfileTwoCosine[5]
<2>_{j,-j} with nested loops.
Double_t ComputeThreeNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
Double_t ComputeFourNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
TProfile * fProfileCosineFourNestedLoops[6]
<4>_{i,j,-i,-j}, ij: kl,km,lm,kn,ln,mn.
TComplex CalculateRecursion(Int_t n, Int_t *harmonic, Int_t mult=1, Int_t skip=0)
TH1D * fHistoPhiPreCuts
Pseudorapidity distribution.
const char Option_t
Definition: External.C:48
TProfile * fProfileCosineFourParticles[6]
<3>_{h,i,j} with nested loops.
bool Bool_t
Definition: External.C:53
TH1D * fHistoEtaPostCuts
Transverse momentum distribution.
TProfile * fProfileCosineThreeNestedLoops[4]
<3>_{h,i,j}, h: (k+l),(k-l) (first two), j: (l-k),(-k-l) (last two).
Bool_t fProcessBothKineAndReco
Type of analysis: MC or AOD.
TProfile * fProfileTwoCosineNestedLoops[5]
<<2>_{i,-i}<2>_{j,-j}> for ij: kl,km,lm, and <<2>_{i,-i}<2>_{j,-j}<2>_{h,-h}>, ijh: klm and <<4>_{k...
virtual void GSCfullAnalysis(long long nParticles, Double_t angles[], Double_t weights[])
Double_t ComputeTwoNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)