AliPhysics  9b6b435 (9b6b435)
AliAnalysisTaskTwoMultiCorrelations.cxx
Go to the documentation of this file.
1 
2 /*************************************************************************
3 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16 
17 //--------------------------------------------------------------------------------------//
18 // Analysis task for the computation of the multiTrack correlations with different //
19 // flow harmonics v_n. This script can takes the Monte Carlo simulations data (e.g. //
20 // HIJING), as well as the experimental Pb-Pb data taken by the ALICE experiment. //
21 // The current script computes the multiTrack correlators using the method of the //
22 // Q-vectors for a maximum of 6 different harmonics and 8 particles). //
23 // //
24 // Author: Cindy Mordasini (cindy.mordasini@cern.ch) //
25 // Version: 12.02.2019 //
26 //--------------------------------------------------------------------------------------//
27 
28 #include "Riostream.h"
29 #include <vector>
31 #include "AliLog.h"
32 #include "AliAODEvent.h"
33 #include "AliAODInputHandler.h"
34 #include "AliMCEvent.h"
35 #include "AliMCEventHandler.h"
36 #include "AliAnalysisManager.h"
37 #include "AliMultSelection.h"
38 #include "AliAODVertex.h"
39 #include "AliMCVertex.h"
40 #include "TList.h"
41 #include "TH1D.h"
42 #include "TFile.h"
43 #include "TComplex.h"
44 #include "TMath.h"
45 
46 using std::cout;
47 using std::endl;
48 
50 
51 //======================================================================================//
52 // Main methods for AliAnalysisTaskSE. //
53 //--------------------------------------------------------------------------------------//
55  AliAnalysisTaskSE(name),
56 // General parameters.
57  fMaxNumberCorrelations(8),
58  fMaxFlowHarmonic(6),
59  fNumberHarmonicsInSC(3),
60  fUseParticleWeights(kFALSE),
61  fComputeNestedLoops(kFALSE),
62 // Structure of the output file.
63  fMainList(NULL),
64  fControlListEventCuts(NULL),
65  fControlListTrackCuts(NULL),
66  fListCorrelations(NULL),
67  fListTwoParticles(NULL),
68  fListThreeParticles(NULL),
69  fListFourParticles(NULL),
70  fListSixParticles(NULL),
71 // Control histograms for the distribution of observables for the event selection.
72  fHistoCentrality(NULL),
73  fHistoNumberOfTracksBeforeCuts(NULL),
74  fHistoNumberOfTracksAfterEventCuts(NULL),
75  fHistoVertexXBeforeCuts(NULL),
76  fHistoVertexYBeforeCuts(NULL),
77  fHistoVertexZBeforeCuts(NULL),
78  fHistoVertexXAfterEventCuts(NULL),
79  fHistoVertexYAfterEventCuts(NULL),
80  fHistoVertexZAfterEventCuts(NULL),
81 // Control histrograms for the distribution of observables for the track selection.
82  fHistoNumberOfTracksAfterAllCuts(NULL),
83  fHistoPtBeforeCuts(NULL),
84  fHistoPtAfterCuts(NULL),
85  fHistoEtaBeforeCuts(NULL),
86  fHistoEtaAfterCuts(NULL),
87  fHistoPhiBeforeCuts(NULL),
88  fHistoPhiAfterCuts(NULL),
89  fHistoTPCClustersBeforeCuts(NULL),
90  fHistoTPCClustersAfterCuts(NULL),
91  fHistoDCAXYBeforeCuts(NULL),
92  fHistoDCAZBeforeCuts(NULL),
93  fHistoDCAXYAfterCuts(NULL),
94  fHistoDCAZAfterCuts(NULL),
95 // TProfiles with the final multiTrack correlations.
96  fProfileCosineEightParticles(NULL),
97 // Type of analysis.
98  fAnalysisType(NULL),
99  fProcessBothKineAndReco(kFALSE),
100  fProcessOnlyKine(kFALSE),
101  fProcessOnlyReco(kFALSE),
102 // Determination of the centrality.
103  fCentralitySelection(NULL),
104  fUseSPDForCentrality(kFALSE),
105  fUseVZeroForCentrality(kFALSE),
106  fCentralityMin(0.),
107  fCentralityMax(100.),
108 // Event selection.
109  fCutOnVertexX(kFALSE),
110  fVertexMinX(-44.),
111  fVertexMaxX(-44.),
112  fCutOnVertexY(kFALSE),
113  fVertexMinY(-44.),
114  fVertexMaxY(-44.),
115  fCutOnVertexZ(kFALSE),
116  fVertexMinZ(-10.),
117  fVertexMaxZ(10.),
118 // Track selection.
119  fPtMin(0.2),
120  fPtMax(5.),
121  fEtaMin(-0.8),
122  fEtaMax(0.8),
123  fNumberOfTPCMin(70),
124  fDCAxyMax(3.2),
125  fDCAzMax(2.4),
126 // Harmonics.
127  fHarmonicOne(2),
128  fHarmonicTwo(-2),
129  fHarmonicThree(3),
130  fHarmonicFour(-3),
131  fHarmonicFive(4),
132  fHarmonicSix(-4),
133  fHarmonicSeven(0),
134  fHarmonicEight(0)
135 {
136 /* Constructor of the class. */
137  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
138 
139 // Create the main list.
140  fMainList = new TList();
141  fMainList->SetName("outputAnalysis");
142  fMainList->SetOwner(kTRUE);
143 
144 // Define the input and output slots.
145  DefineOutput(1, TList::Class());
146 
147 // Initialise the fQvectors to zero.
148  InitialiseArraysOfQvectors();
149 
150 // Initialise the pointers of the TProfiles to NULL.
151  InitialiseArraysOfTProfiles();
152 
153 // Use the particle weights?
154  if(useParticleWeights)
155  {
156  /* TBA, but not needed for LHC10h */
157  }
158 } // End: AliAnalysisTaskTwoMultiCorrelations().
159 
160 //======================================================================================//
163 // General parameters.
164  fMaxNumberCorrelations(8),
165  fMaxFlowHarmonic(6),
166  fNumberHarmonicsInSC(3),
167  fUseParticleWeights(kFALSE),
168  fComputeNestedLoops(kFALSE),
169 // Structure of the output file.
170  fMainList(NULL),
171  fControlListEventCuts(NULL),
172  fControlListTrackCuts(NULL),
173  fListCorrelations(NULL),
174  fListTwoParticles(NULL),
175  fListThreeParticles(NULL),
176  fListFourParticles(NULL),
177  fListSixParticles(NULL),
178 // Control histograms for the distribution of observables for the event selection.
179  fHistoCentrality(NULL),
180  fHistoNumberOfTracksBeforeCuts(NULL),
181  fHistoNumberOfTracksAfterEventCuts(NULL),
182  fHistoVertexXBeforeCuts(NULL),
183  fHistoVertexYBeforeCuts(NULL),
184  fHistoVertexZBeforeCuts(NULL),
185  fHistoVertexXAfterEventCuts(NULL),
186  fHistoVertexYAfterEventCuts(NULL),
187  fHistoVertexZAfterEventCuts(NULL),
188 // Control histrograms for the distribution of observables for the track selection.
189  fHistoNumberOfTracksAfterAllCuts(NULL),
190  fHistoPtBeforeCuts(NULL),
191  fHistoPtAfterCuts(NULL),
192  fHistoEtaBeforeCuts(NULL),
193  fHistoEtaAfterCuts(NULL),
194  fHistoPhiBeforeCuts(NULL),
195  fHistoPhiAfterCuts(NULL),
196  fHistoTPCClustersBeforeCuts(NULL),
197  fHistoTPCClustersAfterCuts(NULL),
198  fHistoDCAXYBeforeCuts(NULL),
199  fHistoDCAZBeforeCuts(NULL),
200  fHistoDCAXYAfterCuts(NULL),
201  fHistoDCAZAfterCuts(NULL),
202 // TProfiles with the final multiTrack correlations.
203  fProfileCosineEightParticles(NULL),
204 // Type of analysis.
205  fAnalysisType(NULL),
206  fProcessBothKineAndReco(kFALSE),
207  fProcessOnlyKine(kFALSE),
208  fProcessOnlyReco(kFALSE),
209 // Determination of the centrality.
210  fCentralitySelection(NULL),
211  fUseSPDForCentrality(kFALSE),
212  fUseVZeroForCentrality(kFALSE),
213  fCentralityMin(0.),
214  fCentralityMax(100.),
215 // Event selection.
216  fCutOnVertexX(kFALSE),
217  fVertexMinX(-44.),
218  fVertexMaxX(-44.),
219  fCutOnVertexY(kFALSE),
220  fVertexMinY(-44.),
221  fVertexMaxY(-44.),
222  fCutOnVertexZ(kFALSE),
223  fVertexMinZ(-10.),
224  fVertexMaxZ(10.),
225 // Track selection.
226  fPtMin(0.2),
227  fPtMax(5.),
228  fEtaMin(-0.8),
229  fEtaMax(0.8),
230  fNumberOfTPCMin(70),
231  fDCAxyMax(3.2),
232  fDCAzMax(2.4),
233 // Harmonics.
234  fHarmonicOne(2),
235  fHarmonicTwo(-2),
236  fHarmonicThree(3),
237  fHarmonicFour(-3),
238  fHarmonicFive(4),
239  fHarmonicSix(-4),
240  fHarmonicSeven(0),
241  fHarmonicEight(0)
242 {
243 /* Dummy constructor of the class. */
244  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
245 
246 // Initialise the fQvectors to zero.
248 
249 // Initialise the pointers of the TProfiles to NULL.
251 
252 } // End: AliAnalysisTaskTwoMultiCorrelations().
253 
254 //======================================================================================//
256 {
257 /* Destructor of the class. */
258  if (fMainList) {delete fMainList;}
259 
260 } // End: ~AliAnalysisTaskTwoMultiCorrelations().
261 
262 //======================================================================================//
264 {
265 /* Steps to do at the beginning of the analysis. */
266 // Avoid name clashes.
267  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
268  TH1::AddDirectory(kFALSE);
269 
270 // Book all the lists.
271  this->BookAllLists();
272 
273 // Book the graphs in all the secondary lists.
274  this->BookControlListEventCuts();
275  this->BookControlListTrackCuts();
276  this->BookListCorrelations();
277 
278 // Initialise the pointers for the TStrings.
279  fAnalysisType = new TString();
281 
282 // Still avoid name clashes.
283  TH1::AddDirectory(oldHistAddStatus);
284  PostData(1,fMainList);
285 
286 } // End: UserCreateOutputObjects().
287 
288 //======================================================================================//
290 {
291 /* Deal with the main analysis for each event. */
292 // TBA: non-unit particle weights for non-uniform acceptance.
293  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::UserExec(Option_t *)";
294 
295 // Select the type of analysis (MC/AOD) from TaskSE.
296  AliMCEvent *currentMCEvent = MCEvent(); // Pointer to a Monte Carlo event.
297  AliAODEvent *currentAODEvent = dynamic_cast<AliAODEvent*>(InputEvent()); // Pointer to an AOD event.
299  {
300  Fatal(sMethodName.Data(), "ERROR: only one fProcess must be kTRUE in 'SetAnalysisType'.");
301  }
302  else if (fProcessBothKineAndReco) {*fAnalysisType = "MC_AOD";}
303  else if (fProcessOnlyKine) {*fAnalysisType = "MC";}
304  else if (fProcessOnlyReco) {*fAnalysisType = "AOD";}
305 
306 // Call the specific analysis.
307  if (fAnalysisType->EqualTo("AOD")) {AODanalysis(currentAODEvent);}
308  else if (fAnalysisType->EqualTo("MC")) {MCanalysis(currentMCEvent);}
309  else if (fAnalysisType->EqualTo("MC_AOD")) {Fatal(sMethodName.Data(),"ERROR: TBA.");}
310  else {Fatal(sMethodName.Data(), "ERROR: fAnalysisType not defined.");}
311 
312 // PostData.
313  PostData(1, fMainList);
314 
315 } // End: UserExec().
316 
317 //======================================================================================//
318 
320 {
321 /* Post-analysis steps done once the run over the events is done. */
322 // Access the main list.
323  fMainList = (TList*)GetOutputData(1);
324  if (!fMainList) {exit(1);}
325 
326 // Create the output file and save the main list inside.
327  TFile *outputFile = new TFile("AnalysisResults.root", "RECREATE");
328  fMainList->Write(fMainList->GetName(),TObject::kSingleKey);
329  delete outputFile;
330 
331 } // End: Terminate().
332 
333 //======================================================================================//
334 // Methods called in the constructor. //
335 //--------------------------------------------------------------------------------------//
337 {
338 /* Initialisation of the Q-vectors to zero. */
339  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
340  {
341  for (Int_t iPower = 0; iPower < 9; iPower++)
342  {
343  fQvectors[iHarmo][iPower] = TComplex(0.,0.);
344  }
345  }
346 
347 } // End: InitialiseArraysOfQvectors().
348 
349 //======================================================================================//
351 {
352 /* Initialisation of the pointers to the TProfiles. */
353  for (Int_t i = 0; i < 6; i++)
354  {
355  fProfileCosineTwoParticles[i] = NULL;
357  fProfileCosineFourParticles[i] = NULL;
359  }
360 
361  for (Int_t j = 0; j < 5; j++)
362  {
363  fProfileTwoCosine[j] = NULL;
365  }
366 
367  for (Int_t k = 0; k < 4; k++)
368  {
371  fProfileCosineSixParticles[k] = NULL;
372  }
373 
374 } // End: InitialiseArraysOfTProfiles().
375 
376 //======================================================================================//
377 // Methods called in UserCreateOutputObjects. //
378 //--------------------------------------------------------------------------------------//
380 {
381 /* Book all the lists. */
382  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::BookAllLists()";
383  if (!fMainList) {Fatal(sMethodName.Data(), "Error: 'fMainList' is NULL.");}
384 
385 // List with the observables used to select the events.
387  fControlListEventCuts->SetName("ControlListEventCuts");
388  fControlListEventCuts->SetOwner(kTRUE);
389  fMainList->Add(fControlListEventCuts);
390 
391 // List with the observables used to select the tracks.
393  fControlListTrackCuts->SetName("ControlListTrackCuts");
394  fControlListTrackCuts->SetOwner(kTRUE);
396 
397 // List with the final multiTrack correlations, separated into sublists according to the number of particles involved in the correlations.
398  fListCorrelations = new TList();
399  fListCorrelations->SetName("ListCorrelations");
400  fListCorrelations->SetOwner(kTRUE);
402 
403  fListTwoParticles = new TList();
404  fListTwoParticles->SetName("ListTwoParticles");
405  fListTwoParticles->SetOwner(kTRUE);
407 
408  fListThreeParticles = new TList();
409  fListThreeParticles->SetName("ListThreeParticles");
410  fListThreeParticles->SetOwner(kTRUE);
412 
413  fListFourParticles = new TList();
414  fListFourParticles->SetName("ListFourParticles");
415  fListFourParticles->SetOwner(kTRUE);
417 
418  fListSixParticles = new TList();
419  fListSixParticles->SetName("ListSixParticles");
420  fListSixParticles->SetOwner(kTRUE);
422 
423 } // End: BookAllLists().
424 
425 //======================================================================================//
427 {
428 /* Book the histograms with the observables for the event selection. */
429 // Centrality distribution.
430  fHistoCentrality = new TH1D("fHistoCentrality", "Distribution of the centrality", 100, 0., 100.);
431  fHistoCentrality->SetStats(kTRUE);
432  fHistoCentrality->GetXaxis()->SetTitle("Centrality percentile");
434 
435 // Multiplicity distribution.
436  fHistoNumberOfTracksBeforeCuts = new TH1D("fHistoNumberOfTracksBeforeCuts", "Distribution of the number of tracks before any selection", 50000, 0., 50000.);
437  fHistoNumberOfTracksBeforeCuts->SetStats(kTRUE);
438  fHistoNumberOfTracksBeforeCuts->GetXaxis()->SetTitle("Number of tracks");
440 
441  fHistoNumberOfTracksAfterEventCuts = new TH1D("fHistoNumberOfTracksAfterEventCuts", "Distribution of the number of tracks after the event selection", 50000, 0., 50000.);
442  fHistoNumberOfTracksAfterEventCuts->SetStats(kTRUE);
443  fHistoNumberOfTracksAfterEventCuts->GetXaxis()->SetTitle("Number of tracks");
445 
446 // x-position of the PV distribution.
447  fHistoVertexXBeforeCuts = new TH1D("fHistoVertexXBeforeCuts", "Distribution of PV_{x} before any selection", 1000, -20., 20.);
448  fHistoVertexXBeforeCuts->SetStats(kTRUE);
449  fHistoVertexXBeforeCuts->GetXaxis()->SetTitle("PV_{x}");
451 
452  fHistoVertexXAfterEventCuts = new TH1D("fHistoVertexXAfterEventCuts", "Distribution of PV_{x} after the event selection", 1000, -20., 20.);
453  fHistoVertexXAfterEventCuts->SetStats(kTRUE);
454  fHistoVertexXAfterEventCuts->GetXaxis()->SetTitle("PV_{x}");
456 
457 // y-position of the PV distribution.
458  fHistoVertexYBeforeCuts = new TH1D("fHistoVertexYBeforeCuts", "Distribution of PV_{y} before any selection", 1000, -20., 20.);
459  fHistoVertexYBeforeCuts->SetStats(kTRUE);
460  fHistoVertexYBeforeCuts->GetXaxis()->SetTitle("PV_{y}");
462 
463  fHistoVertexYAfterEventCuts = new TH1D("fHistoVertexYAfterEventCuts", "Distribution of PV_{y} after the event selection", 1000, -20., 20.);
464  fHistoVertexYAfterEventCuts->SetStats(kTRUE);
465  fHistoVertexYAfterEventCuts->GetXaxis()->SetTitle("PV_{y}");
467 
468 // z-position of the PV distribution.
469  fHistoVertexZBeforeCuts = new TH1D("fHistoVertexZBeforeCuts", "Distribution of PV_{z} before any selection", 1000, -20., 20.);
470  fHistoVertexZBeforeCuts->SetStats(kTRUE);
471  fHistoVertexZBeforeCuts->GetXaxis()->SetTitle("PV_{z}");
473 
474  fHistoVertexZAfterEventCuts = new TH1D("fHistoVertexZAfterEventCuts", "Distribution of PV_{z} after the event selection", 1000, -20., 20.);
475  fHistoVertexZAfterEventCuts->SetStats(kTRUE);
476  fHistoVertexZAfterEventCuts->GetXaxis()->SetTitle("PV_{z}");
478 
479 } // End: BookControlListEventCuts().
480 
481 //======================================================================================//
483 {
484 /* Book the histograms with the observables for the track selection. */
485 // Multiplicity distribution.
486  fHistoNumberOfTracksAfterAllCuts = new TH1D("fHistoNumberOfTracksAfterAllCuts", "Distribution of the number of tracks after the track selection", 50000, 0., 50000.);
487  fHistoNumberOfTracksAfterAllCuts->SetStats(kTRUE);
488  fHistoNumberOfTracksAfterAllCuts->GetXaxis()->SetTitle("Number of tracks");
490 
491 // Transverse momentum distribution.
492  fHistoPtBeforeCuts = new TH1D("fHistoPtBeforeCuts", "Distribution of p_{T} before the track selection", 1000, 0., 20.);
493  fHistoPtBeforeCuts->SetStats(kTRUE);
494  fHistoPtBeforeCuts->GetXaxis()->SetTitle("p_{T}");
496 
497  fHistoPtAfterCuts = new TH1D("fHistoPtAfterCuts", "Distribution of p_{T} after the track selection", 1000, 0., 20.);
498  fHistoPtAfterCuts->SetStats(kTRUE);
499  fHistoPtAfterCuts->GetXaxis()->SetTitle("p_{T}");
501 
502 // Pseudorapidity distribution.
503  fHistoEtaBeforeCuts = new TH1D("fHistoEtaBeforeCuts", "Distribution of #eta before the track selection", 1000, -1., 1.);
504  fHistoEtaBeforeCuts->SetStats(kTRUE);
505  fHistoEtaBeforeCuts->GetXaxis()->SetTitle("#eta");
507 
508  fHistoEtaAfterCuts = new TH1D("fHistoEtaAfterCuts", "Distribution of #eta after the track selection", 1000, -1., 1.);
509  fHistoEtaAfterCuts->SetStats(kTRUE);
510  fHistoEtaAfterCuts->GetXaxis()->SetTitle("#eta");
512 
513 // Azimuthal angles distribution.
514  fHistoPhiBeforeCuts = new TH1D("fHistoPhiBeforeCuts", "Distributiion of #phi before the track selection", 1000, 0., 6.3);
515  fHistoPhiBeforeCuts->SetStats(kTRUE);
516  fHistoPhiBeforeCuts->GetXaxis()->SetTitle("#phi");
518 
519  fHistoPhiAfterCuts = new TH1D("fHistoPhiAfterCuts", "Distribution of #phi after the track selection", 1000, 0., 6.3);
520  fHistoPhiAfterCuts->SetStats(kTRUE);
521  fHistoPhiAfterCuts->GetXaxis()->SetTitle("#phi");
523 
524 // TPC cluster distribution.
525  fHistoTPCClustersBeforeCuts = new TH1D("fHistoTPCClustersBeforeCuts", "Distribution of the number of TPC clusters before the track selection", 1000, 0., 160.);
526  fHistoTPCClustersBeforeCuts->SetStats(kTRUE);
527  fHistoTPCClustersBeforeCuts->GetXaxis()->SetTitle("TPC clusters");
529 
530  fHistoTPCClustersAfterCuts = new TH1D("fHistoTPCClustersAfterCuts", "Distribution of the number of TPC clusters after the track selection", 1000, 0., 160.);
531  fHistoTPCClustersAfterCuts->SetStats(kTRUE);
532  fHistoTPCClustersAfterCuts->GetXaxis()->SetTitle("TPC clusters");
534 
535 // xy-plane of the DCA distribution.
536  fHistoDCAXYBeforeCuts = new TH1D("fHistoDCAXYBeforeCuts", "Distribution of DCA_{xy} before the track selection", 1000, 0., 10.);
537  fHistoDCAXYBeforeCuts->SetStats(kTRUE);
538  fHistoDCAXYBeforeCuts->GetXaxis()->SetTitle("DCA_{xy}");
540 
541  fHistoDCAXYAfterCuts = new TH1D("fHistoDCAXYAfterCuts", "Distribution of DCA_{xy} after the track selection", 1000, 0., 10.);
542  fHistoDCAXYAfterCuts->SetStats(kTRUE);
543  fHistoDCAXYAfterCuts->GetXaxis()->SetTitle("DCA_{xy}");
545 
546 // z-coordinate of the DCA distribution.
547  fHistoDCAZBeforeCuts = new TH1D("fHistoDCAZBeforeCuts", "Distribution of DCA_{z} before the track selection", 1000, 0., 10.);
548  fHistoDCAZBeforeCuts->SetStats(kTRUE);
549  fHistoDCAZBeforeCuts->GetXaxis()->SetTitle("DCA_{z}");
551 
552  fHistoDCAZAfterCuts = new TH1D("fHistoDCAZAfterCuts", "Distribution of DCA_{z} after the track selection", 1000, 0., 10.);
553  fHistoDCAZAfterCuts->SetStats(kTRUE);
554  fHistoDCAZAfterCuts->GetXaxis()->SetTitle("DCA_{z}");
556 
557 } // End: BookControlListTrackCuts().
558 
559 //======================================================================================//
561 {
562 /* Book the TProfiles with the results for the multiTrack correlations. */
563 // For 2-particle correlations.
564  for (Int_t t = 0; t < 6; t++)
565  {
566  fProfileCosineTwoParticles[t] = new TProfile("", "", 1, 0., 1.);
567  fProfileCosineTwoParticles[t]->SetStats(kTRUE);
568  fProfileCosineTwoParticles[t]->Sumw2();
569  fProfileCosineTwoParticles[t]->SetName(Form("fProfileCosineTwoParticles_%d", t));
571 
573  {
574  fProfileCosineTwoNestedLoops[t] = new TProfile("", "", 1, 0., 1.);
575  fProfileCosineTwoNestedLoops[t]->SetStats(kTRUE);
576  fProfileCosineTwoNestedLoops[t]->Sumw2();
577  fProfileCosineTwoNestedLoops[t]->SetName(Form("fProfileCosineTwoNestedLoops_%d", t));
579  }
580  }
581 
582  for (Int_t c = 0; c < 5; c++)
583  {
584  fProfileTwoCosine[c] = new TProfile("", "", 1, 0., 1.);
585  fProfileTwoCosine[c]->SetStats(kTRUE);
586  fProfileTwoCosine[c]->Sumw2();
587  fProfileTwoCosine[c]->SetName(Form("fProfileTwoCosine_%d", c));
589 
591  {
592  fProfileTwoCosineNestedLoops[c] = new TProfile("", "", 1, 0., 1.);
593  fProfileTwoCosineNestedLoops[c]->SetStats(kTRUE);
595  fProfileTwoCosineNestedLoops[c]->SetName(Form("fProfileTwoCosineNestedLoops_%d", c));
597  }
598  }
599 
600 // For 3-particle correlations.
601  for (Int_t th = 0; th < 4; th++)
602  {
603  fProfileCosineThreeParticles[th] = new TProfile("","", 1, 0., 1.);
604  fProfileCosineThreeParticles[th]->SetStats(kTRUE);
605  fProfileCosineThreeParticles[th]->Sumw2();
606  fProfileCosineThreeParticles[th]->SetName(Form("fProfileCosineThreeParticles_%d", th));
608 
610  {
611  fProfileCosineThreeNestedLoops[th] = new TProfile("","", 1, 0., 1.);
612  fProfileCosineThreeNestedLoops[th]->SetStats(kTRUE);
613  fProfileCosineThreeNestedLoops[th]->Sumw2();
614  fProfileCosineThreeNestedLoops[th]->SetName(Form("fProfileCosineThreeNestedLoops_%d", th));
616  }
617  }
618 
619 // For 4-particle correlations.
620  for (Int_t f = 0; f < 6; f++)
621  {
622  fProfileCosineFourParticles[f] = new TProfile("", "", 1, 0., 1.);
623  fProfileCosineFourParticles[f]->SetStats(kTRUE);
624  fProfileCosineFourParticles[f]->Sumw2();
625  fProfileCosineFourParticles[f]->SetName(Form("fProfileCosineFourParticles_%d", f));
627 
629  {
630  fProfileCosineFourNestedLoops[f] = new TProfile("", "", 1, 0., 1.);
631  fProfileCosineFourNestedLoops[f]->SetStats(kTRUE);
632  fProfileCosineFourNestedLoops[f]->Sumw2();
633  fProfileCosineFourNestedLoops[f]->SetName(Form("fProfileCosineFourNestedLoops_%d", f));
635  }
636  }
637 
638 // For 6- and 8-particle correlations.
639  for (Int_t s = 0; s < 4; s++)
640  {
641  fProfileCosineSixParticles[s] = new TProfile("", "", 1, 0., 1.);
642  fProfileCosineSixParticles[s]->SetStats(kTRUE);
643  fProfileCosineSixParticles[s]->Sumw2();
644  fProfileCosineSixParticles[s]->SetName(Form("fProfileCosineSixParticles_%d", s));
646  }
647 
648  fProfileCosineEightParticles = new TProfile("", "", 1, 0., 1.);
649  fProfileCosineEightParticles->SetStats(kTRUE);
651  fProfileCosineEightParticles->SetName("fProfileCosineEightParticles");
653 
654 } // End: BookListCorrelations().
655 
656 //======================================================================================//
657 // Methods called in UserExec(Option_t *).
658 //--------------------------------------------------------------------------------------//
660 {
661 /* Compute the multiTrack correlations for an AOD event. */
662  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::AODanalysis(AliAODEvent *aAODevent)";
663 
664 // Determine the detector to use for the centrality estimations.
666  {
667  Fatal(sMethodName.Data(), "ERROR: only one detector must be selected in 'SetCentralityEstimation'.");
668  }
669  else if (fUseSPDForCentrality) {*fCentralitySelection = "CL1";}
670  else if (fUseVZeroForCentrality) {*fCentralitySelection = "V0M";}
671 
672 // Determine if the event belongs to this centrality range.
673  AliMultSelection *ams = (AliMultSelection*)aAODevent->FindListObject("MultSelection");
674  if (!ams) {return;} // Protection against NULL pointer.
675  if (ams->GetMultiplicityPercentile(Form("%s", fCentralitySelection->Data())) >= fCentralityMin && ams->GetMultiplicityPercentile(Form("%s", fCentralitySelection->Data())) < fCentralityMax)
676  {
677  fHistoCentrality->Fill(ams->GetMultiplicityPercentile(Form("%s", fCentralitySelection->Data())));
678  }
679  else {return;} // This event does not belong to this centrality range.
680 
681 // Fill the histogram for the distribution of the number of tracks before selection.
682  long long numberOfTracksBeforeCuts = aAODevent->GetNumberOfTracks(); // Number of tracks before any selection.
683  fHistoNumberOfTracksBeforeCuts->Fill(numberOfTracksBeforeCuts);
684 
685 //......................................................................................//
686 /* Application of the event selection. */
687 // Cuts on the PV position.
688  AliAODVertex *avtx = (AliAODVertex*)aAODevent->GetPrimaryVertex(); // Position of the PV.
689  fHistoVertexXBeforeCuts->Fill(avtx->GetX());
690  fHistoVertexYBeforeCuts->Fill(avtx->GetY());
691  fHistoVertexZBeforeCuts->Fill(avtx->GetZ());
692 
693  if (fCutOnVertexX)
694  {
695  if ((avtx->GetX() < fVertexMinX) || (avtx->GetX() > fVertexMaxX)) {return;}
696  }
697  if (fCutOnVertexY)
698  {
699  if ((avtx->GetY() < fVertexMinY) || (avtx->GetY() > fVertexMaxY)) {return;}
700  }
701  if (fCutOnVertexZ)
702  {
703  if ((avtx->GetZ() < fVertexMinZ) || (avtx->GetZ() > fVertexMaxZ)) {return;}
704  }
705 
706 // Fill the histograms after the event selection.
707  fHistoVertexXAfterEventCuts->Fill(avtx->GetX());
708  fHistoVertexYAfterEventCuts->Fill(avtx->GetY());
709  fHistoVertexZAfterEventCuts->Fill(avtx->GetZ());
710 
711 //......................................................................................//
712 /* Application of the track selection. */
713  long long numberOfTracksAfterEventCuts = aAODevent->GetNumberOfTracks(); // Number of tracks after the event selection and before the track selection.
714  long long nParticlesAfterCuts = 0; // Number of tracks after the track selection (= multiplicity).
715  Int_t *hasTrackPassedTheCuts = new Int_t[numberOfTracksAfterEventCuts](); // Flag if a track passed the selection (1) or not (0).
716 
717  fHistoNumberOfTracksAfterEventCuts->Fill(numberOfTracksAfterEventCuts);
718 
719  for (long long iTrack = 0; iTrack < numberOfTracksAfterEventCuts; iTrack++)
720  {
721  // Get a pointer to the current track.
722  AliAODTrack *currentTrack = dynamic_cast<AliAODTrack*>(aAODevent->GetTrack(iTrack));
723  if (!currentTrack) {continue;} // Protection against NULL pointer.
724  if (!currentTrack->TestFilterBit(128)){continue;} // Filter bit 128 denotes TPC-only tracks.
725 
726  // Get all the observables for the selection.
727  Double_t ptBeforeCuts = currentTrack->Pt(); // Transverse momentum.
728  Double_t etaBeforeCuts = currentTrack->Eta(); // Pseudorapidity.
729  Double_t phiBeforeCuts = currentTrack->Phi(); // Azimuthal angle.
730  Int_t numberOfTPCclustersBeforeCuts = currentTrack->GetTPCncls(); // Number of TPC clusters.
731  Double_t DCAxBeforeCuts = currentTrack->XAtDCA(); // x-coordinate of the DCA.
732  Double_t DCAyBeforeCuts = currentTrack->YAtDCA(); // y-coordinate of the DCA.
733  Double_t DCAzBeforeCuts = currentTrack->ZAtDCA(); // z-coordinate of the DCA.
734 
735  Double_t DCAxyBeforeCuts = TMath::Sqrt((DCAxBeforeCuts*DCAxBeforeCuts) + (DCAyBeforeCuts*DCAyBeforeCuts));
736 
737  // Fill the control histograms of some observables before the cuts.
738  fHistoPtBeforeCuts->Fill(ptBeforeCuts);
739  fHistoEtaBeforeCuts->Fill(etaBeforeCuts);
740  fHistoPhiBeforeCuts->Fill(phiBeforeCuts);
741  fHistoTPCClustersBeforeCuts->Fill(numberOfTPCclustersBeforeCuts);
742  fHistoDCAXYBeforeCuts->Fill(DCAxyBeforeCuts);
743  fHistoDCAZBeforeCuts->Fill(DCAzBeforeCuts);
744 
745  // Apply the track selection to the current track.
746  if (CreateTrackSelection(ptBeforeCuts, etaBeforeCuts, numberOfTPCclustersBeforeCuts, DCAxyBeforeCuts, DCAzBeforeCuts))
747  {
748  nParticlesAfterCuts++;
749  hasTrackPassedTheCuts[iTrack] = 1;
750  }
751  else {hasTrackPassedTheCuts[iTrack] = 0;}
752  }
753 
754  fHistoNumberOfTracksAfterAllCuts->Fill(nParticlesAfterCuts);
755 
756 // Define the observables after the track selection for the AOD analysis.
757  Double_t *phi = new Double_t[nParticlesAfterCuts](); // Azimuthal angle.
758  Double_t *particleWeights = new Double_t[nParticlesAfterCuts](); // Particle weight.
759  Int_t indexNewArray = 0; // New index of the selected particle after the selection.
760  Double_t pt = 0.; // Transverse momentum.
761  Double_t eta = 0.; // Pseudorapidity.
762  Int_t numberOfTPCclusters = 0; // Number of TPC clusters.
763  Double_t DCAx = 0.; // x-coordinate of the DCA.
764  Double_t DCAy = 0.; // y-coordinate of the DCA.
765  Double_t DCAz = 0.; // z-coordinate of the DCA.
766  Double_t DCAxy = 0.; // xy-plane of the DCA.
767 
768  if (nParticlesAfterCuts <= 2*fNumberHarmonicsInSC) {return;} // Do the analysis only if the number of particles after the cuts is at least equal to the maximum number of particles needed for the studied GSC.
769 
770  for (long long iiTrack = 0; iiTrack < numberOfTracksAfterEventCuts; iiTrack++)
771  {
772  // Pointer to the current particle.
773  AliAODTrack *keptParticle = dynamic_cast<AliAODTrack*>(aAODevent->GetTrack(iiTrack));
774  if (!keptParticle) {continue;} // Protection against NULL pointer.
775  if(!keptParticle->TestFilterBit(128)) {continue;} // Filter bit 128 denotes TPC-only tracks.
776 
777  if (hasTrackPassedTheCuts[iiTrack] == 1) // The particle passed the selection.
778  {
779  // Get the values of the observables for the current particle.
780  phi[indexNewArray] = keptParticle->Phi();
781  pt = keptParticle->Pt();
782  eta = keptParticle->Eta();
783  numberOfTPCclusters = keptParticle->GetTPCncls();
784  DCAx = keptParticle->XAtDCA();
785  DCAy = keptParticle->YAtDCA();
786  DCAz = keptParticle->ZAtDCA();
787 
788  DCAxy = TMath::Sqrt((DCAx*DCAx) + (DCAy*DCAy));
789 
790  particleWeights[indexNewArray] = 1.;
791  //if (fUseParticleWeights) {continue;} // TBA
792  //else {particleWeights[indexNewArray] = 1.;}
793 
794  // Fill the control histograms after the track selection.
795  fHistoPhiAfterCuts->Fill(phi[indexNewArray]);
796  fHistoPtAfterCuts->Fill(pt);
797  fHistoEtaAfterCuts->Fill(eta);
798  fHistoTPCClustersAfterCuts->Fill(numberOfTPCclusters);
799  fHistoDCAXYAfterCuts->Fill(DCAxy);
800  fHistoDCAZAfterCuts->Fill(DCAz);
801 
802  indexNewArray++;
803  }
804  else {continue;}
805  }
806 
807 // Calculate all the Q-vectors for the current event.
808  CalculateQvectors(nParticlesAfterCuts, phi, particleWeights);
809 
810 // Compute all the multiTrack correlations for the current event.
811  GSCfullAnalysis(nParticlesAfterCuts, phi, particleWeights);
812 
813 // Reset everything to zero for the next event.
814  delete [] phi;
815  delete [] particleWeights;
816  indexNewArray = 0;
817  pt = 0.;
818  eta = 0.;
819  numberOfTPCclusters = 0;
820  DCAx = 0.;
821  DCAy = 0.;
822  DCAz = 0.;
823  DCAxy = 0.;
824 
825 } // End: AODanalysis().
826 
827 //======================================================================================//
829 {
830 /* Compute the multiTrack correlations for a MC event. */
831  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::MCanalysis(AliMCEvent *aMCevent)";
832 
833 // Determine the detector to use for the centrality estimations.
835  {
836  Fatal(sMethodName.Data(), "ERROR: only one detector must be selected in 'SetCentralityEstimation'.");
837  }
838  else if (fUseSPDForCentrality) {*fCentralitySelection = "CL1";}
839  else if (fUseVZeroForCentrality) {*fCentralitySelection = "V0M";}
840 
841 // Determine if the event belongs to this centrality range (for reconstructed particles only).
842  if (!fProcessOnlyKine)
843  {
844  AliMultSelection *ams = (AliMultSelection*)aMCevent->FindListObject("MultSelection");
845  if (!ams) {return;} // Protection against NULL pointer.
846  if (ams->GetMultiplicityPercentile(Form("%s", fCentralitySelection->Data())) >= fCentralityMin && ams->GetMultiplicityPercentile(Form("%s", fCentralitySelection->Data())) < fCentralityMax)
847  {
848  fHistoCentrality->Fill(ams->GetMultiplicityPercentile(Form("%s", fCentralitySelection->Data())));
849  }
850  else {return;} // This event does not belong to this centrality range.
851  }
852 
853 // Fill the histogram for the distribution of the number of tracks before selection.
854  long long numberOfTracksBeforeCuts = aMCevent->GetNumberOfTracks(); // Number of tracks before any selection.
855  fHistoNumberOfTracksBeforeCuts->Fill(numberOfTracksBeforeCuts);
856 
857 //......................................................................................//
858 /* Application of the event selection. */
859 // Cuts on the PV position.
860  AliMCVertex *avtx = (AliMCVertex*)aMCevent->GetPrimaryVertex(); // Position of the PV.
861  fHistoVertexXBeforeCuts->Fill(avtx->GetX());
862  fHistoVertexYBeforeCuts->Fill(avtx->GetY());
863  fHistoVertexZBeforeCuts->Fill(avtx->GetZ());
864 
865  if (fCutOnVertexX)
866  {
867  if ((avtx->GetX() < fVertexMinX) || (avtx->GetX() > fVertexMaxX)) {return;}
868  }
869  if (fCutOnVertexY)
870  {
871  if ((avtx->GetY() < fVertexMinY) || (avtx->GetY() > fVertexMaxY)) {return;}
872  }
873  if (fCutOnVertexZ)
874  {
875  if ((avtx->GetZ() < fVertexMinZ) || (avtx->GetZ() > fVertexMaxZ)) {return;}
876  }
877 
878 // Fill the histograms after the event selection.
879  fHistoVertexXAfterEventCuts->Fill(avtx->GetX());
880  fHistoVertexYAfterEventCuts->Fill(avtx->GetY());
881  fHistoVertexZAfterEventCuts->Fill(avtx->GetZ());
882 
883 //......................................................................................//
884 /* Application of the track selection. */
885  long long numberOfTracksAfterEventCuts = aMCevent->GetNumberOfTracks(); // Number of tracks after the event selection and before the track selection.
886  long long nParticlesAfterCuts = 0; // Number of tracks after the track selection (= multiplicity).
887  Int_t *hasTrackPassedTheCuts = new Int_t[numberOfTracksAfterEventCuts](); // Flag if a track passed the selection (1) or not (0).
888 
889  fHistoNumberOfTracksAfterEventCuts->Fill(numberOfTracksAfterEventCuts);
890 
891  for (long long iTrack = 0; iTrack < numberOfTracksAfterEventCuts; iTrack++)
892  {
893  // Get a pointer to the current track.
894  AliAODMCParticle *currentTrack = dynamic_cast<AliAODMCParticle*>(aMCevent->GetTrack(iTrack));
895  if (!currentTrack) {continue;} // Protection against NULL pointer.
896 
897  // Get all the observables for the selection.
898  Double_t ptBeforeCuts = currentTrack->Pt(); // Transverse momentum.
899  Double_t etaBeforeCuts = currentTrack->Eta(); // Pseudorapidity.
900  Double_t phiBeforeCuts = currentTrack->Phi(); // Azimuthal angle.
901 
902  // Fill the control histograms of some observables before the cuts.
903  fHistoPtBeforeCuts->Fill(ptBeforeCuts);
904  fHistoEtaBeforeCuts->Fill(etaBeforeCuts);
905  fHistoPhiBeforeCuts->Fill(phiBeforeCuts);
906 
907  // Apply the track selection to the current track.
908  if ((fPtMin <= ptBeforeCuts) && (ptBeforeCuts <= fPtMax) && (fEtaMin <= etaBeforeCuts) && (etaBeforeCuts <= fEtaMax))
909  {
910  nParticlesAfterCuts++;
911  hasTrackPassedTheCuts[iTrack] = 1;
912  }
913  else {hasTrackPassedTheCuts[iTrack] = 0;}
914  }
915 
916  fHistoNumberOfTracksAfterAllCuts->Fill(nParticlesAfterCuts);
917 
918 // Define the observables after the track selection for the AOD analysis.
919  Double_t *phi = new Double_t[nParticlesAfterCuts](); // Azimuthal angle.
920  Double_t *particleWeights = new Double_t[nParticlesAfterCuts](); // Particle weight.
921  Int_t indexNewArray = 0; // New index of the selected particle after the selection.
922  Double_t pt = 0.; // Transverse momentum.
923  Double_t eta = 0.; // Pseudorapidity.
924 
925  if (nParticlesAfterCuts <= 2*fNumberHarmonicsInSC) {return;} // Do the analysis only if the number of particles after the cuts is at least equal to the maximum number of particles needed for the studied GSC.
926 
927  for (long long iiTrack = 0; iiTrack < numberOfTracksAfterEventCuts; iiTrack++)
928  {
929  // Pointer to the current particle.
930  AliAODMCParticle *keptParticle = dynamic_cast<AliAODMCParticle*>(aMCevent->GetTrack(iiTrack));
931  if (!keptParticle) {continue;} // Protection against NULL pointer.
932 
933  if (hasTrackPassedTheCuts[iiTrack] == 1) // The particle passed the selection.
934  {
935  // Get the values of the observables for the current particle.
936  phi[indexNewArray] = keptParticle->Phi();
937  pt = keptParticle->Pt();
938  eta = keptParticle->Eta();
939 
940  particleWeights[indexNewArray] = 1.;
941  //if (fUseParticleWeights) {continue;} // TBA
942  //else {particleWeights[indexNewArray] = 1.;}
943 
944  // Fill the control histograms after the track selection.
945  fHistoPhiAfterCuts->Fill(phi[indexNewArray]);
946  fHistoPtAfterCuts->Fill(pt);
947  fHistoEtaAfterCuts->Fill(eta);
948 
949  indexNewArray++;
950  }
951  else {continue;}
952  }
953 
954 // Calculate all the Q-vectors for the current event.
955  CalculateQvectors(nParticlesAfterCuts, phi, particleWeights);
956 
957 // Compute all the multiTrack correlations for the current event.
958  GSCfullAnalysis(nParticlesAfterCuts, phi, particleWeights);
959 
960 // Reset everything to zero for the next event.
961  delete [] phi;
962  delete [] particleWeights;
963  indexNewArray = 0;
964  pt = 0.;
965  eta = 0.;
966 
967 } // End: MCDanalysis().
968 
969 //======================================================================================//
970 Bool_t AliAnalysisTaskTwoMultiCorrelations::CreateTrackSelection(Double_t currentPt, Double_t currentEta, Int_t currentNumberOfTPC, Double_t currentDCAXY, Double_t currentDCAZ)
971 {
972 /* List the cuts to apply to the particles in a given event. */
973  Bool_t cutOnPt = (fPtMin <= currentPt) && (currentPt <= fPtMax);
974  Bool_t cutOnEta = (fEtaMin <= currentEta) && (currentEta <= fEtaMax);
975  Bool_t cutOnNumberOfTPC = (fNumberOfTPCMin < currentNumberOfTPC);
976  Bool_t cutOnDCAxy = (currentDCAXY < fDCAxyMax);
977  Bool_t cutOnDCAz = (currentDCAZ < fDCAzMax);
978 
979  return cutOnPt && cutOnEta && cutOnNumberOfTPC && cutOnDCAxy && cutOnDCAz;
980 
981 } // End: Bool_t CreateTrackSelection().
982 
983 //======================================================================================//
984 void AliAnalysisTaskTwoMultiCorrelations::CalculateQvectors(long long nParticles, Double_t angles[], Double_t weights[])
985 {
986 /* Calculate all the Q-vectors for a given set of azimuthal angles and particles weights. */
987  Double_t iAngle = 0.; // Angle of the current particle.
988  Double_t iWeight = 1.; // Particle weight of the current particle.
989  Double_t weightToPowerP = 1.; // Particle weight rised to the power p.
990 
991 // Ensure all the Q-vectors are initially zero.
992  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
993  {
994  for (Int_t iPower = 0; iPower < 9; iPower++)
995  {
996  fQvectors[iHarmo][iPower] = TComplex(0.,0.);
997  }
998  }
999 
1000 // Compute the Q-vectors.
1001  for (long long iTrack = 0; iTrack < nParticles; iTrack++)
1002  {
1003  iAngle = angles[iTrack];
1004  if (fUseParticleWeights) {iWeight = weights[iTrack];}
1005  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
1006  {
1007  for (Int_t iPower = 0; iPower < 9; iPower++)
1008  {
1009  if (fUseParticleWeights) {weightToPowerP = TMath::Power(iWeight, iPower);}
1010  fQvectors[iHarmo][iPower] += TComplex(weightToPowerP*TMath::Cos(iHarmo*iAngle), weightToPowerP*TMath::Sin(iHarmo*iAngle));
1011  }
1012  }
1013  }
1014 } // End: CalculateQvectors().
1015 
1016 //======================================================================================//
1017 void AliAnalysisTaskTwoMultiCorrelations::GSCfullAnalysis(long long nParticles, Double_t angles[], Double_t weights[])
1018 {
1019 /* Compute all the needed multiTrack correlators (with possible cross-check from nested loops) for the chosen harmonics. */
1020  Int_t harmonicsArray[8] = {fHarmonicOne, fHarmonicTwo, fHarmonicThree, fHarmonicFour, fHarmonicFive, fHarmonicSix, fHarmonicSeven, fHarmonicEight}; // Harmonics (n_1,... n_8).
1021 
1022 // Compute the denominator of each case of correlators.
1023  Int_t twoZerosArray[2] = {0}; // For 2p correlations.
1024  Double_t twoParticleDenominator = CalculateRecursion(2, twoZerosArray).Re();
1025  Int_t threeZerosArray[3] = {0}; // For 3p correlations.
1026  Double_t threeParticleDenominator = CalculateRecursion(3, threeZerosArray).Re();
1027  Int_t fourZerosArray[4] = {0}; // For 4p correlations.
1028  Double_t fourParticleDenominator = CalculateRecursion(4, fourZerosArray).Re();
1029  Int_t sixZerosArray[6] = {0}; // For 6p correlations.
1030  Double_t sixParticleDenominator = CalculateRecursion(6, sixZerosArray).Re();
1031  Int_t eightZerosArray[8] = {0}; // For 8p correlations.
1032  Double_t eightParticleDenominator = CalculateRecursion(8, eightZerosArray).Re();
1033  eightParticleDenominator++; // Dummy line since this variable is not useful yet.
1034 
1035 //**************************************************************************************//
1036 // Compute the 2-particle correlations.
1037  Int_t firstHarmonicArray[2] = {harmonicsArray[0], harmonicsArray[1]}; // Array with (k,-k).
1038  Int_t secondHarmonicArray[2] = {harmonicsArray[2], harmonicsArray[3]}; // Array with (l,-l).
1039  Int_t thirdHarmonicArray[2] = {harmonicsArray[4], harmonicsArray[5]}; // Array with (m,-m).
1040  Int_t fourthHarmonicArray[2] = {harmonicsArray[6], harmonicsArray[7]}; // Array with (n,-n).
1041 
1042  Double_t tempoThirdHarmoNested = 0.;
1043 
1045  TComplex firstHarmonicComplex = (CalculateRecursion(2, firstHarmonicArray))/twoParticleDenominator; // Complex value of <2>_{k,-k}.
1046  Double_t firstHarmonicValue = firstHarmonicComplex.Re(); // Value of <2>_{k,-k}.
1047  fProfileCosineTwoParticles[0]->Fill(0.5, firstHarmonicValue, twoParticleDenominator);
1048  fProfileCosineTwoParticles[0]->SetTitle(Form("#LT2#GT_{k,-k}, k = %d", firstHarmonicArray[0]));
1049  fProfileCosineTwoParticles[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", firstHarmonicArray[0], firstHarmonicArray[1]));
1050 
1052  TComplex secondHarmonicComplex = (CalculateRecursion(2, secondHarmonicArray))/twoParticleDenominator; // Complex value of <2>_{l,-l}.
1053  Double_t secondHarmonicValue = secondHarmonicComplex.Re(); // Value of <2>_{l,-l}.
1054  fProfileCosineTwoParticles[1]->Fill(0.5, secondHarmonicValue, twoParticleDenominator);
1055  fProfileCosineTwoParticles[1]->SetTitle(Form("#LT2#GT_{l,-l}, l = %d", secondHarmonicArray[0]));
1056  fProfileCosineTwoParticles[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", secondHarmonicArray[0], secondHarmonicArray[1]));
1057 
1059  TComplex thirdHarmonicComplex = TComplex(0.,0.); // Complex value of <2>_{m,-m}.
1060  Double_t thirdHarmonicValue = 0.; // Value of <2>_{m,-m}.
1061  if (fNumberHarmonicsInSC >= 3)
1062  {
1063  thirdHarmonicComplex = (CalculateRecursion(2, thirdHarmonicArray))/twoParticleDenominator;
1064  thirdHarmonicValue = thirdHarmonicComplex.Re();
1065  fProfileCosineTwoParticles[2]->Fill(0.5, thirdHarmonicValue, twoParticleDenominator);
1066  fProfileCosineTwoParticles[2]->SetTitle(Form("#LT2#GT_{m,-m}, m = %d", thirdHarmonicArray[0]));
1067  fProfileCosineTwoParticles[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", thirdHarmonicArray[0], thirdHarmonicArray[1]));
1068  } // End: if (fNumberHarmonicsInSC >= 3)
1069 
1071  TComplex fourthHarmonicComplex = TComplex(0.,0.); // Complex value of <2>_{n,-n}.
1072  Double_t fourthHarmonicValue = 0.; // Value of <2>_{n,-n}.
1073  if (fNumberHarmonicsInSC == 4)
1074  {
1075  fourthHarmonicComplex = (CalculateRecursion(2, fourthHarmonicArray))/twoParticleDenominator;
1076  fourthHarmonicValue = fourthHarmonicComplex.Re();
1077  fProfileCosineTwoParticles[3]->Fill(0.5, fourthHarmonicValue, twoParticleDenominator);
1078  fProfileCosineTwoParticles[3]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", fourthHarmonicArray[0]));
1079  fProfileCosineTwoParticles[3]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", fourthHarmonicArray[0], fourthHarmonicArray[1]));
1080  } // End: if (fNumberHarmonicsInSC == 4)
1081 
1083  Int_t sumTwoHarmonicsArray[2] = {firstHarmonicArray[0]+secondHarmonicArray[0], firstHarmonicArray[1]+secondHarmonicArray[1]}; // Array with (k+l,-k-l).
1084  TComplex sumTwoHarmonicsComplex = TComplex(0.,0.); // Complex value of <2>_{k+l,-k-l}.
1085  Double_t sumTwoHarmonicsValue = 0.; // <2>_{k+l,-k-l}.
1086  if (fNumberHarmonicsInSC == 2)
1087  {
1088  sumTwoHarmonicsComplex = (CalculateRecursion(2, sumTwoHarmonicsArray))/twoParticleDenominator;
1089  sumTwoHarmonicsValue = sumTwoHarmonicsComplex.Re();
1090  fProfileCosineTwoParticles[4]->Fill(0.5, sumTwoHarmonicsValue, twoParticleDenominator);
1091  fProfileCosineTwoParticles[4]->SetTitle(Form("#LT2#GT_{k+l,-k-l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
1092  fProfileCosineTwoParticles[4]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", sumTwoHarmonicsArray[0], sumTwoHarmonicsArray[1]));
1093  } // End: if (fNumberHarmonicsInSC == 2)
1094 
1096  Int_t diffTwoHarmonicsArray[2] = {firstHarmonicArray[0]+secondHarmonicArray[1], firstHarmonicArray[1]+secondHarmonicArray[0]}; // Array with (k-l,-k+l).
1097  TComplex diffTwoHarmonicsComplex = TComplex(0.,0.); // Complex value of <2>_{k-l,-k+l}.
1098  Double_t diffTwoHarmonicsValue = 0.; // <2>_{k-l,-k+l}.
1099  if (fNumberHarmonicsInSC == 2)
1100  {
1101  diffTwoHarmonicsComplex = (CalculateRecursion(2, diffTwoHarmonicsArray))/twoParticleDenominator;
1102  diffTwoHarmonicsValue = diffTwoHarmonicsComplex.Re();
1103  fProfileCosineTwoParticles[5]->Fill(0.5, diffTwoHarmonicsValue, twoParticleDenominator);
1104  fProfileCosineTwoParticles[5]->SetTitle(Form("#LT2#GT_{k-l,-k+l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
1105  fProfileCosineTwoParticles[5]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d,%d}", diffTwoHarmonicsArray[0], diffTwoHarmonicsArray[1]));
1106  } // End: if (fNumberHarmonicsInSC == 2)
1107 
1109  Double_t doubleCosineValueFirstCase = firstHarmonicValue*secondHarmonicValue;
1110  fProfileTwoCosine[0]->Fill(0.5, doubleCosineValueFirstCase);
1111  fProfileTwoCosine[0]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{l,-l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
1112  fProfileTwoCosine[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], secondHarmonicArray[0], secondHarmonicArray[1]));
1113 
1115  Double_t doubleCosineValueSecondCase = 0.;
1116  if (fNumberHarmonicsInSC >= 3)
1117  {
1118  doubleCosineValueSecondCase = firstHarmonicValue*thirdHarmonicValue;
1119  fProfileTwoCosine[1]->Fill(0.5, doubleCosineValueSecondCase);
1120  fProfileTwoCosine[1]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{m,-m}, k = %d, m = %d", firstHarmonicArray[0], thirdHarmonicArray[0]));
1121  fProfileTwoCosine[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
1122  } // End: if (fNumberHarmonicsInSC >= 3)
1123 
1125  Double_t doubleCosineValueThirdCase = 0.;
1126  if (fNumberHarmonicsInSC >= 3)
1127  {
1128  doubleCosineValueThirdCase = secondHarmonicValue*thirdHarmonicValue;
1129  fProfileTwoCosine[2]->Fill(0.5, doubleCosineValueThirdCase);
1130  fProfileTwoCosine[2]->SetTitle(Form("#LT2#GT_{l,-l}#LT2#GT_{m,-m}, l = %d, m = %d", secondHarmonicArray[0], thirdHarmonicArray[0]));
1131  fProfileTwoCosine[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", secondHarmonicArray[0], secondHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
1132  } // End: if (fNumberHarmonicsInSC >= 3)
1133 
1135  Double_t doubleCosineValueFourthCase = 0.;
1136  if (fNumberHarmonicsInSC >= 3)
1137  {
1138  doubleCosineValueFourthCase = firstHarmonicValue*secondHarmonicValue*thirdHarmonicValue;
1139  fProfileTwoCosine[3]->Fill(0.5, doubleCosineValueFourthCase);
1140  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]));
1141  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]));
1142  } // End: if (fNumberHarmonicsInSC >= 3)
1143 
1145  if (fComputeNestedLoops)
1146  {
1147  // 1st harmonic: <2>_{k,-k}.
1148  Double_t firstHarmonicValueNested = ComputeTwoNestedLoops(nParticles, firstHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[0]);
1149  fProfileCosineTwoNestedLoops[0]->SetTitle(Form("#LT2#GT_{k,-k}, k = %d", firstHarmonicArray[0]));
1150  fProfileCosineTwoNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", firstHarmonicArray[0], firstHarmonicArray[1]));
1151 
1152  // 2nd harmonic: <2>_{l,-l}.
1153  Double_t secondHarmonicValueNested = ComputeTwoNestedLoops(nParticles, secondHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[1]);
1154  fProfileCosineTwoNestedLoops[1]->SetTitle(Form("#LT2#GT_{l,-l}, l = %d", secondHarmonicArray[0]));
1155  fProfileCosineTwoNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", secondHarmonicArray[0], secondHarmonicArray[1]));
1156 
1157  // 3rd harmonic: <2>_{m,-m}.
1158  Double_t thirdHarmonicValueNested = 0.;
1159  if (fNumberHarmonicsInSC >= 3)
1160  {
1161  thirdHarmonicValueNested = ComputeTwoNestedLoops(nParticles, thirdHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[2]);
1162  fProfileCosineTwoNestedLoops[2]->SetTitle(Form("#LT2#GT_{m,-m}, m = %d", thirdHarmonicArray[0]));
1163  fProfileCosineTwoNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", thirdHarmonicArray[0], thirdHarmonicArray[1]));
1164  } // End: if (fNumberHarmonicsInSC >= 3)
1165  tempoThirdHarmoNested = thirdHarmonicValueNested;
1166 
1167  // 4th harmonic: <2>_{n,-n}.
1168  Double_t fourthHarmonicValueNested = 0.;
1169  fourthHarmonicValueNested++; // Dummy line since this variable is not useful yet.
1170  if (fNumberHarmonicsInSC == 4)
1171  {
1172  fourthHarmonicValueNested = ComputeTwoNestedLoops(nParticles, fourthHarmonicArray, angles, weights, fProfileCosineTwoNestedLoops[3]);
1173  fProfileCosineTwoNestedLoops[3]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", fourthHarmonicArray[0]));
1174  fProfileCosineTwoNestedLoops[3]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", fourthHarmonicArray[0], fourthHarmonicArray[1]));
1175  } // End: if (fNumberHarmonicsInSC == 4)
1176 
1177  // Sum first two harmonics: <2>_{k+l,-k-l}.
1178  Double_t sumTwoHarmonicsValueNested = 0.;
1179  sumTwoHarmonicsValueNested++; // Dummy line since this variable is not useful yet.
1180  if (fNumberHarmonicsInSC == 2)
1181  {
1182  sumTwoHarmonicsValueNested = ComputeTwoNestedLoops(nParticles, sumTwoHarmonicsArray, angles, weights, fProfileCosineTwoNestedLoops[4]);
1183  fProfileCosineTwoNestedLoops[4]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", sumTwoHarmonicsArray[0]));
1184  fProfileCosineTwoNestedLoops[4]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", sumTwoHarmonicsArray[0], sumTwoHarmonicsArray[1]));
1185  } // End: if (fNumberHarmonicsInSC == 2)
1186 
1187  // Difference first two harmonics: <2>_{k-l,-k+l}.
1188  Double_t diffTwoHarmonicsValueNested = 0.;
1189  diffTwoHarmonicsValueNested++; // Dummy line since this variable is not useful yet.
1190  if (fNumberHarmonicsInSC == 2)
1191  {
1192  diffTwoHarmonicsValueNested = ComputeTwoNestedLoops(nParticles, diffTwoHarmonicsArray, angles, weights, fProfileCosineTwoNestedLoops[5]);
1193  fProfileCosineTwoNestedLoops[5]->SetTitle(Form("#LT2#GT_{n,-n}, n = %d", diffTwoHarmonicsArray[0]));
1194  fProfileCosineTwoNestedLoops[5]->GetYaxis()->SetTitle(Form("#LT#LT2#GT#GT_{%d, %d}", diffTwoHarmonicsArray[0], diffTwoHarmonicsArray[1]));
1195  } // End: if (fNumberHarmonicsInSC == 2)
1196 
1197  // <<2>_{k,-k}<2>_{l,-l}>.
1198  Double_t doubleCosineFirstCaseNested = firstHarmonicValueNested*secondHarmonicValueNested;
1199  fProfileTwoCosineNestedLoops[0]->Fill(0.5, doubleCosineFirstCaseNested);
1200  fProfileTwoCosineNestedLoops[0]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{l,-l}, k = %d, l = %d", firstHarmonicArray[0], secondHarmonicArray[0]));
1201  fProfileTwoCosineNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], secondHarmonicArray[0], secondHarmonicArray[1]));
1202 
1203  // <2>_{k,-k}<2>_{m,-m}.
1204  Double_t doubleCosineSecondCaseNested = 0.;
1205  if (fNumberHarmonicsInSC >= 3)
1206  {
1207  doubleCosineSecondCaseNested = firstHarmonicValueNested*thirdHarmonicValueNested;
1208  fProfileTwoCosineNestedLoops[1]->Fill(0.5, doubleCosineSecondCaseNested);
1209  fProfileTwoCosineNestedLoops[1]->SetTitle(Form("#LT2#GT_{k,-k}#LT2#GT_{m,-m}, k = %d, m = %d", firstHarmonicArray[0], thirdHarmonicArray[0]));
1210  fProfileTwoCosineNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", firstHarmonicArray[0], firstHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
1211  } // End: if (fNumberHarmonicsInSC >= 3)
1212 
1213  // <2>_{l,-l}<2>_{m,-m}.
1214  Double_t doubleCosineThirdCaseNested = 0.;
1215  if (fNumberHarmonicsInSC >= 3)
1216  {
1217  doubleCosineThirdCaseNested = secondHarmonicValueNested*thirdHarmonicValueNested;
1218  fProfileTwoCosineNestedLoops[2]->Fill(0.5, doubleCosineThirdCaseNested);
1219  fProfileTwoCosineNestedLoops[2]->SetTitle(Form("#LT2#GT_{l,-l}#LT2#GT_{m,-m}, l = %d, m = %d", secondHarmonicArray[0], thirdHarmonicArray[0]));
1220  fProfileTwoCosineNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT2#GT_{%d,%d}#LT2#GT_{%d,%d}#GT", secondHarmonicArray[0], secondHarmonicArray[1], thirdHarmonicArray[0], thirdHarmonicArray[1]));
1221  } // End: if (fNumberHarmonicsInSC >= 3)
1222 
1223  // <2>_{k,-k}<2>_{l,-l}<2>_{m,-m}.
1224  Double_t doubleCosineFourthCaseNested = 0.;
1225  if (fNumberHarmonicsInSC >= 3)
1226  {
1227  doubleCosineFourthCaseNested = firstHarmonicValueNested*secondHarmonicValueNested*thirdHarmonicValueNested;
1228  fProfileTwoCosineNestedLoops[3]->Fill(0.5, doubleCosineFourthCaseNested);
1229  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]));
1230  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]));
1231  } // End: if (fNumberHarmonicsInSC >= 3)
1232 
1233  // Reset of the observables for the two nested loops.
1234  firstHarmonicValueNested = 0.;
1235  secondHarmonicValueNested = 0.;
1236  thirdHarmonicValueNested = 0.;
1237  fourthHarmonicValueNested = 0.;
1238  sumTwoHarmonicsValueNested = 0.;
1239  diffTwoHarmonicsValueNested = 0.;
1240  doubleCosineFirstCaseNested = 0.;
1241  doubleCosineSecondCaseNested = 0.;
1242  doubleCosineThirdCaseNested = 0.;
1243  doubleCosineFourthCaseNested = 0.;
1244  } // End: if (fComputeNestedLoops)
1245 
1246 //**************************************************************************************//
1247 // Compute the 3-particle correlations.
1249  Int_t threeParticleCaseOneArray[3] = {harmonicsArray[0]+harmonicsArray[2], harmonicsArray[1], harmonicsArray[3]};
1250  TComplex threeParticleCaseOneComplex = TComplex(0.,0.); // Complex value of <3>_{k+l,-k,-l}.
1251  Double_t threeParticleCaseOneValue = 0.; // Value of <3>_{k+l,-k,-l}.
1252  if (fNumberHarmonicsInSC == 2)
1253  {
1254  threeParticleCaseOneComplex = (CalculateRecursion(3, threeParticleCaseOneArray))/threeParticleDenominator;
1255  threeParticleCaseOneValue = threeParticleCaseOneComplex.Re();
1256  fProfileCosineThreeParticles[0]->Fill(0.5, threeParticleCaseOneValue, threeParticleDenominator);
1257  fProfileCosineThreeParticles[0]->SetTitle(Form("#LT3#GT_{k+l,-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1258  fProfileCosineThreeParticles[0]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseOneArray[0], threeParticleCaseOneArray[1], threeParticleCaseOneArray[2]));
1259  } // End: if (fNumberHarmonicsInSC == 2)
1260 
1262  Int_t threeParticleCaseTwoArray[3] = {harmonicsArray[0]+harmonicsArray[3], harmonicsArray[1], harmonicsArray[2]};
1263  TComplex threeParticleCaseTwoComplex = TComplex(0.,0.); // Complex value of <3>_{k+l,-k,-l}.
1264  Double_t threeParticleCaseTwoValue = 0.; // Value of <3>_{k-l,-k,+l}.
1265  if (fNumberHarmonicsInSC == 2)
1266  {
1267  threeParticleCaseTwoComplex = (CalculateRecursion(3, threeParticleCaseTwoArray))/threeParticleDenominator;
1268  threeParticleCaseTwoValue = threeParticleCaseTwoComplex.Re();
1269  fProfileCosineThreeParticles[1]->Fill(0.5, threeParticleCaseTwoValue, threeParticleDenominator);
1270  fProfileCosineThreeParticles[1]->SetTitle(Form("#LT3#GT_{k-l,-k,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1271  fProfileCosineThreeParticles[1]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseTwoArray[0], threeParticleCaseTwoArray[1], threeParticleCaseTwoArray[2]));
1272  } // End: if (fNumberHarmonicsInSC == 2)
1273 
1275  Int_t threeParticleCaseThreeArray[3] = {harmonicsArray[0], harmonicsArray[2]+harmonicsArray[1], harmonicsArray[3]};
1276  TComplex threeParticleCaseThreeComplex = TComplex(0.,0.); // Complex value of <3>_{k,l-k,-l}.
1277  Double_t threeParticleCaseThreeValue = 0.; // Value of <3>_{k,l-k,-l}.
1278  if (fNumberHarmonicsInSC == 2)
1279  {
1280  threeParticleCaseThreeComplex = (CalculateRecursion(3, threeParticleCaseThreeArray))/threeParticleDenominator;
1281  threeParticleCaseThreeValue = threeParticleCaseThreeComplex.Re();
1282  fProfileCosineThreeParticles[2]->Fill(0.5, threeParticleCaseThreeValue, threeParticleDenominator);
1283  fProfileCosineThreeParticles[2]->SetTitle(Form("#LT3#GT_{k,l-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1284  fProfileCosineThreeParticles[2]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseThreeArray[0], threeParticleCaseThreeArray[1], threeParticleCaseThreeArray[2]));
1285  } // End: if (fNumberHarmonicsInSC == 2)
1286 
1288  Int_t threeParticleCaseFourArray[3] = {harmonicsArray[0], harmonicsArray[1]+harmonicsArray[3], harmonicsArray[2]};
1289  TComplex threeParticleCaseFourComplex = TComplex(0.,0.); // Complex value of <3>_{k,-k-l,l}.
1290  Double_t threeParticleCaseFourValue = 0.; // Value of <3>_{k,-k-l,l}.
1291  if (fNumberHarmonicsInSC == 2)
1292  {
1293  threeParticleCaseFourComplex = (CalculateRecursion(3, threeParticleCaseFourArray))/threeParticleDenominator;
1294  threeParticleCaseFourValue = threeParticleCaseFourComplex.Re();
1295  fProfileCosineThreeParticles[3]->Fill(0.5, threeParticleCaseFourValue, threeParticleDenominator);
1296  fProfileCosineThreeParticles[3]->SetTitle(Form("#LT3#GT_{k,-k-l,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1297  fProfileCosineThreeParticles[3]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseFourArray[0], threeParticleCaseFourArray[1], threeParticleCaseFourArray[2]));
1298  } // End: if (fNumberHarmonicsInSC == 2)
1299 
1301  if (fComputeNestedLoops)
1302  {
1303  // Case one: <3>_{k+l,-k,-l}.
1304  Double_t threeParticleCaseOneNested = 0.;
1305  threeParticleCaseOneNested++; // Dummy line since this variable is not useful yet.
1306  if (fNumberHarmonicsInSC == 2)
1307  {
1308  threeParticleCaseOneNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseOneArray, angles, weights, fProfileCosineThreeNestedLoops[0]);
1309  fProfileCosineThreeNestedLoops[0]->SetTitle(Form("#LT3#GT_{k+l,-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1310  fProfileCosineThreeNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseOneArray[0], threeParticleCaseOneArray[1], threeParticleCaseOneArray[2]));
1311  } // End: if (fNumberHarmonicsInSC == 2)
1312 
1313  // Case two: <3>_{k-l,-k,+l}.
1314  Double_t threeParticleCaseTwoNested = 0.;
1315  threeParticleCaseTwoNested++; // Dummy line since this variable is not useful yet.
1316  if (fNumberHarmonicsInSC == 2)
1317  {
1318  threeParticleCaseTwoNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseTwoArray, angles, weights, fProfileCosineThreeNestedLoops[1]);
1319  fProfileCosineThreeNestedLoops[1]->SetTitle(Form("#LT3#GT_{k-l,-k,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1320  fProfileCosineThreeNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseTwoArray[0], threeParticleCaseTwoArray[1], threeParticleCaseTwoArray[2]));
1321  } // End: if (fNumberHarmonicsInSC == 2)
1322 
1323  // Case three: <3>_{k,l-k,-l}.
1324  Double_t threeParticleCaseThreeNested = 0.;
1325  threeParticleCaseThreeNested++; // Dummy line since this variable is not useful yet.
1326  if (fNumberHarmonicsInSC == 2)
1327  {
1328  threeParticleCaseThreeNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseThreeArray, angles, weights, fProfileCosineThreeNestedLoops[2]);
1329  fProfileCosineThreeNestedLoops[2]->SetTitle(Form("#LT3#GT_{k,l-k,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1330  fProfileCosineThreeNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseThreeArray[0], threeParticleCaseThreeArray[1], threeParticleCaseThreeArray[2]));
1331  } // End: if (fNumberHarmonicsInSC == 2)
1332 
1333  // Case four: <3>_{k,-k-l,l}.
1334  Double_t threeParticleCaseFourNested = 0.;
1335  threeParticleCaseFourNested++; // Dummy line since this variable is not useful yet.
1336  if (fNumberHarmonicsInSC == 2)
1337  {
1338  threeParticleCaseFourNested = ComputeThreeNestedLoops(nParticles, threeParticleCaseFourArray, angles, weights, fProfileCosineThreeNestedLoops[3]);
1339  fProfileCosineThreeNestedLoops[3]->SetTitle(Form("#LT3#GT_{k,-k-l,l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1340  fProfileCosineThreeNestedLoops[3]->GetYaxis()->SetTitle(Form("#LT#LT3#GT#GT_{%d,%d,%d}", threeParticleCaseFourArray[0], threeParticleCaseFourArray[1], threeParticleCaseFourArray[2]));
1341  } // End: if (fNumberHarmonicsInSC == 2)
1342 
1343  // Reset of the observables for the three nested loops.
1344  threeParticleCaseOneNested = 0.;
1345  threeParticleCaseTwoNested = 0.;
1346  threeParticleCaseThreeNested = 0.;
1347  threeParticleCaseFourNested = 0.;
1348  } // End: if (fComputeNestedLoops)
1349 
1350 //**************************************************************************************//
1351 // Compute the 4-particle correlations.
1353  Int_t fourParticleCaseOneArray[4] = {harmonicsArray[0], harmonicsArray[1], harmonicsArray[2], harmonicsArray[3]};
1354  TComplex fourParticleCaseOneComplex = (CalculateRecursion(4, fourParticleCaseOneArray))/fourParticleDenominator; // Complex value of <4>_{k,-k,l,-l}.
1355  Double_t fourParticleCaseOneValue = fourParticleCaseOneComplex.Re(); // Value of <4>_{k,-k,l,-l}.
1356  fProfileCosineFourParticles[0]->Fill(0.5, fourParticleCaseOneValue, fourParticleDenominator);
1357  fProfileCosineFourParticles[0]->SetTitle(Form("#LT4#GT_{k,-k,l,-l}, k = %d, l = %d", fourParticleCaseOneArray[0], fourParticleCaseOneArray[2]));
1358  fProfileCosineFourParticles[0]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseOneArray[0], fourParticleCaseOneArray[1],fourParticleCaseOneArray[2],fourParticleCaseOneArray[3]));
1359 
1361  Int_t fourParticleCaseTwoArray[4] = {harmonicsArray[0], harmonicsArray[1], harmonicsArray[4], harmonicsArray[5]};
1362  TComplex fourParticleCaseTwoComplex = TComplex(0.,0.); // Complex value of <4>_{k,-k,m,-m}.
1363  Double_t fourParticleCaseTwoValue = 0.; // Value of <4>_{k,-k,m,-m}.
1364  if (fNumberHarmonicsInSC >= 3)
1365  {
1366  fourParticleCaseTwoComplex = (CalculateRecursion(4, fourParticleCaseTwoArray))/fourParticleDenominator;
1367  fourParticleCaseTwoValue = fourParticleCaseTwoComplex.Re();
1368  fProfileCosineFourParticles[1]->Fill(0.5, fourParticleCaseTwoValue, fourParticleDenominator);
1369  fProfileCosineFourParticles[1]->SetTitle(Form("#LT4#GT_{k,-k,m,-m}, k = %d, m = %d", fourParticleCaseTwoArray[0], fourParticleCaseTwoArray[2]));
1370  fProfileCosineFourParticles[1]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseTwoArray[0], fourParticleCaseTwoArray[1],fourParticleCaseTwoArray[2],fourParticleCaseTwoArray[3]));
1371  } // End: if (fNumberHarmonicsInSC >= 3)
1372 
1374  Int_t fourParticleCaseThreeArray[4] = {harmonicsArray[2], harmonicsArray[3], harmonicsArray[4], harmonicsArray[5]};
1375  TComplex fourParticleCaseThreeComplex = TComplex(0.,0.); // Complex value of <4>_{l,-l,m,-m}.
1376  Double_t fourParticleCaseThreeValue = 0.; // Value of <4>_{l,-l,m,-m}.
1377  if (fNumberHarmonicsInSC >= 3)
1378  {
1379  fourParticleCaseThreeComplex = (CalculateRecursion(4, fourParticleCaseThreeArray))/fourParticleDenominator;
1380  fourParticleCaseThreeValue = fourParticleCaseThreeComplex.Re();
1381  fProfileCosineFourParticles[2]->Fill(0.5, fourParticleCaseThreeValue, fourParticleDenominator);
1382  fProfileCosineFourParticles[2]->SetTitle(Form("#LT4#GT_{l,-l,m,-m}, l = %d, m = %d", fourParticleCaseThreeArray[0], fourParticleCaseThreeArray[2]));
1383  fProfileCosineFourParticles[2]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseThreeArray[0], fourParticleCaseThreeArray[1],fourParticleCaseThreeArray[2],fourParticleCaseThreeArray[3]));
1384  } // End: if (fNumberHarmonicsInSC >= 3)
1385 
1387 
1389  Double_t doubleCosineValueFifthCase = 0.;
1390  if (fNumberHarmonicsInSC >= 3)
1391  {
1392  doubleCosineValueFifthCase = fourParticleCaseOneValue*thirdHarmonicValue;
1393  fProfileTwoCosine[4]->Fill(0.5, doubleCosineValueFifthCase);
1394  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]));
1395  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]));
1396  } // End: if (fNumberHarmonicsInSC >= 3)
1397 
1399  if (fComputeNestedLoops)
1400  {
1401  // Case one: <4>_{k,-k,l,-l}.
1402  Double_t fourParticleCaseOneNested = ComputeFourNestedLoops(nParticles, fourParticleCaseOneArray, angles, weights, fProfileCosineFourNestedLoops[0]);
1403  fProfileCosineFourNestedLoops[0]->SetTitle(Form("#LT4#GT_{k,-k,l,-l}, k = %d, l = %d", harmonicsArray[0], harmonicsArray[2]));
1404  fProfileCosineFourNestedLoops[0]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseOneArray[0], fourParticleCaseOneArray[1], fourParticleCaseOneArray[2],fourParticleCaseOneArray[3]));
1405 
1406  // Case two: <4>_{k,m,-k,-m}.
1407  Double_t fourParticleCaseTwoNested = 0.;
1408  fourParticleCaseTwoNested++; // Dummy line since this variable is not useful yet.
1409  if (fNumberHarmonicsInSC >= 3)
1410  {
1411  fourParticleCaseTwoNested = ComputeFourNestedLoops(nParticles, fourParticleCaseTwoArray, angles, weights, fProfileCosineFourNestedLoops[1]);
1412  fProfileCosineFourNestedLoops[1]->SetTitle(Form("#LT4#GT_{k,-k,m,-m}, k = %d, m = %d", harmonicsArray[0], harmonicsArray[4]));
1413  fProfileCosineFourNestedLoops[1]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseTwoArray[0], fourParticleCaseTwoArray[1], fourParticleCaseTwoArray[2],fourParticleCaseTwoArray[3]));
1414  } // End: if (fNumberHarmonicsInSC >= 3)
1415 
1416  // Case three: <4>_{l,m,-l,-m}.
1417  Double_t fourParticleCaseThreeNested = 0.;
1418  fourParticleCaseThreeNested++; // Dummy line since this variable is not useful yet.
1419  if (fNumberHarmonicsInSC >= 3)
1420  {
1421  fourParticleCaseThreeNested = ComputeFourNestedLoops(nParticles, fourParticleCaseThreeArray, angles, weights, fProfileCosineFourNestedLoops[2]);
1422  fProfileCosineFourNestedLoops[2]->SetTitle(Form("#LT4#GT_{l,-l,m,-m}, l = %d, m = %d", harmonicsArray[2], harmonicsArray[4]));
1423  fProfileCosineFourNestedLoops[2]->GetYaxis()->SetTitle(Form("#LT#LT4#GT#GT_{%d,%d,%d,%d}", fourParticleCaseThreeArray[0], fourParticleCaseThreeArray[1], fourParticleCaseThreeArray[2],fourParticleCaseThreeArray[3]));
1424  } // End: if (fNumberHarmonicsInSC >= 3)
1425 
1426  // TBI: cases with four harmonics...
1427 
1428  // <4>_{k,-k,l,-l}<2>_{m,-m}.
1429  Double_t doubleCosineFifthCaseNested = 0.;
1430  doubleCosineFifthCaseNested++; // Dummy line since this variable is not useful yet.
1431  if (fNumberHarmonicsInSC >= 3)
1432  {
1433  doubleCosineFifthCaseNested = fourParticleCaseOneNested*tempoThirdHarmoNested;
1434  fProfileTwoCosineNestedLoops[4]->Fill(0.5, doubleCosineFifthCaseNested);
1435  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]));
1436  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]));
1437  } // End: if (fNumberHarmonicsInSC >= 3)
1438 
1439  // Reset of the observables for the four nested loops.
1440  fourParticleCaseOneNested = 0.;
1441  fourParticleCaseTwoNested = 0.;
1442  fourParticleCaseThreeNested = 0.;
1443  //
1444  doubleCosineFifthCaseNested = 0.;
1445  } // End: if (fComputeNestedLoops)
1446 
1447 //**************************************************************************************//
1448 // Compute the 6-particle correlations.
1450  Int_t sixParticleCaseOneArray[6] = {harmonicsArray[0], harmonicsArray[1], harmonicsArray[2], harmonicsArray[3], harmonicsArray[4], harmonicsArray[5]};
1451  TComplex sixParticleCaseOneComplex = TComplex(0.,0.); // Complex value of <6>_{k,-k,l,-l,m,-m}.
1452  Double_t sixParticleCaseOneValue = 0.; // Value of <6>_{k,-k,l,-l,m,-m}.
1453  if (fNumberHarmonicsInSC >= 3)
1454  {
1455  sixParticleCaseOneComplex = (CalculateRecursion(6, sixParticleCaseOneArray))/sixParticleDenominator;
1456  sixParticleCaseOneValue = sixParticleCaseOneComplex.Re();
1457  fProfileCosineSixParticles[0]->Fill(0.5, sixParticleCaseOneValue, sixParticleDenominator);
1458  fProfileCosineSixParticles[0]->SetTitle(Form("#LT6#GT_{k,-k,l,-l,m,-m}, k = %d, l = %d, m = %d", sixParticleCaseOneArray[0], sixParticleCaseOneArray[2], sixParticleCaseOneArray[4]));
1459  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]));
1460  } // End: if (fNumberHarmonicsInSC >= 3)
1461 
1462 //**************************************************************************************//
1463 // Compute the 8-particle correlations.
1465 
1466 //**************************************************************************************//
1467 // Reset of the variables to zero.
1468  twoParticleDenominator = 0.;
1469  threeParticleDenominator = 0.;
1470  fourParticleDenominator = 0.;
1471  sixParticleDenominator = 0.;
1472  eightParticleDenominator = 0.;
1473  tempoThirdHarmoNested = 0.;
1474 
1475  firstHarmonicComplex = TComplex(0.,0.);
1476  firstHarmonicValue = 0.;
1477  secondHarmonicComplex = TComplex(0.,0.);
1478  secondHarmonicValue = 0.;
1479  thirdHarmonicComplex = TComplex(0.,0.);
1480  thirdHarmonicValue = 0.;
1481  fourthHarmonicComplex = TComplex(0.,0.);
1482  fourthHarmonicValue = 0.;
1483 
1484  doubleCosineValueFirstCase = 0.,
1485  doubleCosineValueSecondCase = 0.;
1486  doubleCosineValueThirdCase = 0.;
1487 
1488  threeParticleCaseOneComplex = TComplex(0.,0.);
1489  threeParticleCaseOneValue = 0.;
1490  threeParticleCaseTwoComplex = TComplex(0.,0.);
1491  threeParticleCaseTwoValue = 0.;
1492  threeParticleCaseThreeComplex = TComplex(0.,0.);
1493  threeParticleCaseThreeValue = 0.;
1494  threeParticleCaseFourComplex = TComplex(0.,0.);
1495  threeParticleCaseFourValue = 0.;
1496 
1497  fourParticleCaseOneComplex = TComplex(0.,0.);
1498  fourParticleCaseOneValue = 0.;
1499  fourParticleCaseTwoComplex = TComplex(0.,0.);
1500  fourParticleCaseTwoValue = 0.;
1501  fourParticleCaseThreeComplex = TComplex(0.,0.);
1502  fourParticleCaseThreeValue = 0.;
1503  // TBI: cases with four possible harmonics: kn, ln, mn...
1504  doubleCosineValueFifthCase = 0.;
1505 
1506  sixParticleCaseOneComplex = TComplex(0.,0.);
1507  sixParticleCaseOneValue = 0.;
1508 } // End: GSCfullAnalysis().
1509 
1510 //======================================================================================//
1512 {
1513 /* Simplify the use of the Q-vectors. */
1514  if (n >= 0) {return fQvectors[n][p];}
1515  return TComplex::Conjugate(fQvectors[-n][p]);
1516 } // End: Q().
1517 
1518 //======================================================================================//
1520 {
1521 /* Calculate multi-particle correlators by using recursion (an improved faster version) originally developed by Kristjan Gulbrandsen (gulbrand@nbi.dk). */
1522  Int_t nm1 = n-1;
1523  TComplex c(Q(harmonic[nm1], mult));
1524  if (nm1 == 0) return c;
1525  c *= CalculateRecursion(nm1, harmonic);
1526  if (nm1 == skip) return c;
1527 
1528  Int_t multp1 = mult+1;
1529  Int_t nm2 = n-2;
1530  Int_t counter1 = 0;
1531  Int_t hhold = harmonic[counter1];
1532  harmonic[counter1] = harmonic[nm2];
1533  harmonic[nm2] = hhold + harmonic[nm1];
1534  TComplex c2(CalculateRecursion(nm1, harmonic, multp1, nm2));
1535  Int_t counter2 = n-3;
1536 
1537  while (counter2 >= skip) {
1538  harmonic[nm2] = harmonic[counter1];
1539  harmonic[counter1] = hhold;
1540  ++counter1;
1541  hhold = harmonic[counter1];
1542  harmonic[counter1] = harmonic[nm2];
1543  harmonic[nm2] = hhold + harmonic[nm1];
1544  c2 += CalculateRecursion(nm1, harmonic, multp1, counter2);
1545  --counter2;
1546  }
1547  harmonic[nm2] = harmonic[counter1];
1548  harmonic[counter1] = hhold;
1549 
1550  if (mult == 1) return c-c2;
1551  return c-Double_t(mult)*c2;
1552 } // End: CalculateRecursion().
1553 
1554 //======================================================================================//
1555 Double_t AliAnalysisTaskTwoMultiCorrelations::ComputeTwoNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
1556 {
1557 /* Calculate the two-particle correlations using two nested loops. */
1558  Double_t twoParticleCosine = 0.; // cos(k(phi_1 - phi_2)).
1559  Double_t twoParticleWeight = 1.; // Total particle weights.
1560  TProfile twoProfile; // Returns <cos(k(phi1-phi2))>.
1561  twoProfile.Sumw2();
1562 
1563  for (long long m = 0; m < nParticles; m++)
1564  {
1565  for (long long n = 0; n < nParticles; n++)
1566  {
1567  // Remove the autocorrelations m==n.
1568  if (m == n) {continue;}
1569  // Compute cos(k(phi_1-phi_2)).
1570  twoParticleCosine = TMath::Cos(harmonic[0]*angles[m] + harmonic[1]*angles[n]);
1571  twoParticleWeight = weights[m] + weights[n];
1572  profile->Fill(0.5, twoParticleCosine, twoParticleWeight);
1573  twoProfile.Fill(0.5, twoParticleCosine, twoParticleWeight);
1574 
1575  // Reset the variables to default.
1576  twoParticleCosine = 0.;
1577  twoParticleWeight = 1.;
1578  } // End: for (long long n = 0; n < nParticles; n++)
1579  } // End: for (long long m = 0; m < nParticles; m++)
1580 
1581 // Return <cos(k(phi_1-phi_2))>.
1582  return twoProfile.GetBinContent(1);
1583 } // End: ComputeTwoNestedLoops(long long, Int_t*, Double_t[], Double_t[], TProfile*)
1584 
1585 //======================================================================================//
1586 Double_t AliAnalysisTaskTwoMultiCorrelations::ComputeThreeNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
1587 {
1588 /* Calculate the three-particle correlations using three nested loops. */
1589  Double_t threeParticleCosine = 0.; // cos(kphi_1 +lphi_2 +mphi_3).
1590  Double_t threeParticleWeight = 1.; // Total particle weights.
1591  TProfile threeProfile; // Returns <cos(kphi_1 +lphi_2 +mphi_3)>.
1592  threeProfile.Sumw2();
1593 
1594  for (long long k = 0; k < nParticles; k++)
1595  {
1596  for (long long l = 0; l < nParticles; l++)
1597  {
1598  if (k == l) {continue;}
1599 
1600  for (long long m = 0; m < nParticles; m++)
1601  {
1602  if ((k == m) || (l == m)) {continue;}
1603 
1604  // Compute cos(kphi_1 +lphi_2 +mphi_3).
1605  threeParticleCosine = TMath::Cos(harmonic[0]*angles[k] + harmonic[1]*angles[l] + harmonic[2]*angles[m]);
1606  threeParticleWeight = weights[k] + weights[l] + weights[m];
1607  profile->Fill(0.5, threeParticleCosine, threeParticleWeight);
1608  threeProfile.Fill(0.5, threeParticleCosine, threeParticleWeight);
1609 
1610  // Reset the variables to default.
1611  threeParticleCosine = 0.;
1612  threeParticleWeight = 1.;
1613  } // End: for (long long m = 0; m < nParticles; m++)
1614  } // End: for (long long l = 0; l < nParticles; l++)
1615 } // End: for (long long k = 0; k < nParticles; k++)
1616 
1617  // Return <cos(kphi_1 +lphi_2 +mphi_3)>.
1618  return threeProfile.GetBinContent(1);
1619 } // End: ComputeThreeNestedLoops().
1620 
1621 //======================================================================================//
1622 Double_t AliAnalysisTaskTwoMultiCorrelations::ComputeFourNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
1623 {
1624 // Calculate the four-particle correlations using four nested loops.
1625  Double_t fourParticleCosine = 0.; // cos(kphi_1 +lphi_2 +mphi_3 +nphi_4).
1626  Double_t fourParticleWeight = 1.; // Total particle weights.
1627  TProfile fourProfile; // Returns <cos(kphi_1 +lphi_2 +mphi_3 +nphi_4)>.
1628  fourProfile.Sumw2();
1629 
1630  for (long long k = 0; k < nParticles; k++)
1631  {
1632  for (long long l = 0; l < nParticles; l++)
1633  {
1634  if (k == l) {continue;}
1635  for (long long m = 0; m < nParticles; m++)
1636  {
1637  if ((k == m) || (l == m)) {continue;}
1638  for (long long n = 0; n < nParticles; n++)
1639  {
1640  if ((k == n) || (l == n) || (m == n)) {continue;}
1641  // Computate cos(k(phi_1 - phi_2) + l(phi_3-phi_4)).
1642  fourParticleCosine = TMath::Cos(harmonic[0]*angles[k] + harmonic[1]*angles[l] + harmonic[2]*angles[m] + harmonic[3]*angles[n]);
1643  fourParticleWeight = weights[k] + weights[l] + weights[m] + weights[n];
1644  profile->Fill(0.5, fourParticleCosine, fourParticleWeight);
1645  fourProfile.Fill(0.5, fourParticleCosine, fourParticleWeight);
1646 
1647  // Reset the variables to default.
1648  fourParticleCosine = 0.;
1649  fourParticleWeight = 1.;
1650  } // End: for (long long n = 0; n < nParticles; n++)
1651  } // End: for (long long m = 0; m < nParticles; m++)
1652  } // End: for (long long l = 0; l < nParticles; l++)
1653  } // End: for (long long k = 0; k < nParticles; k++)
1654 
1655  // Return <cos(kphi_1 +lphi_2 +mphi_3 +nphi_4)>.
1656  return fourProfile.GetBinContent(1);
1657 } // End: Double_t ComputeFourNestedLoops().
1658 
1659 //======================================================================================//
1660 // Methods called in Terminate(Option_t *).
1661 //--------------------------------------------------------------------------------------//
1662 
1663 //======================================================================================//
1664 
TH1D * fHistoVertexYAfterEventCuts
x-position of the PV after the event selection.
TProfile * fProfileCosineTwoNestedLoops[6]
<2>_{j,-j}, j: k,l,m,n,(k+l),(k-l).
TProfile * fProfileCosineTwoParticles[6]
z-coordinate of the DCA after the track selection.
double Double_t
Definition: External.C:58
Bool_t CreateTrackSelection(Double_t currentPt, Double_t currentEta, Int_t currentNumberOfTPC, Double_t currentDCAXY, Double_t currentDCAZ)
TH1D * fHistoVertexZAfterEventCuts
y-position of the PV after the event selection.
TProfile * fProfileCosineThreeParticles[4]
N<<2>_{i,-i}<2>_{j,-j}> with nested loops.
TH1D * fHistoDCAXYBeforeCuts
Number of TPC clusters after the track selection.
TH1D * fHistoPtBeforeCuts
Number of tracks remaining after both the event and the track selection.
TProfile * fProfileCosineSixParticles[4]
<4>_{i,j,-i,-j} with nested loops.
TH1D * fHistoEtaAfterCuts
Pseudorapidity distribution before the track selection.
TCanvas * c
Definition: TestFitELoss.C:172
TH1D * fHistoVertexXBeforeCuts
Number of tracks remaining after the event selection.
virtual void CalculateQvectors(long long nParticles, Double_t angles[], Double_t weights[])
TH1D * fHistoPhiAfterCuts
Azimuthal angles distribution before the track selection.
TH1D * fHistoTPCClustersBeforeCuts
Azimuthal angles distribution after the track selection.
TH1D * fHistoVertexZBeforeCuts
y-position of the PV before the event selection.
TH1D * fHistoDCAXYAfterCuts
z-coordinate of the DCA before the track selection.
int Int_t
Definition: External.C:63
TH1D * fHistoTPCClustersAfterCuts
Number of TPC clusters before the track selection.
TH1D * fHistoPtAfterCuts
Transverse momentum before the track selection.
TH1D * fHistoEtaBeforeCuts
Transverse momentum distribution after the track selection.
TProfile * fProfileCosineEightParticles
<6>_{h,i,j,-h,-i,-j}, hij: klm,kln,kmn,lmn.
Definition: External.C:212
Bool_t fUseSPDForCentrality
Detector for the centrality estimation: SPD (CL1) or V0 (V0M).
TH1D * fHistoNumberOfTracksBeforeCuts
Centrality distribution.
TString * fAnalysisType
<8>_{k,l,m,n,-k,-l,-m,-n}.
TProfile * fProfileTwoCosine[5]
<2>_{j,-j} with nested loops.
TH1D * fHistoVertexXAfterEventCuts
z-position of the PV before the event selection.
Double_t ComputeThreeNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
TH1D * fHistoPhiBeforeCuts
Pseudorapidity distribution after the track selection.
Double_t ComputeFourNestedLoops(long long nParticles, Int_t *harmonic, Double_t angles[], Double_t weights[], TProfile *profile)
TH1D * fHistoNumberOfTracksAfterAllCuts
z-position of the PV after the event selection.
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)
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 * fHistoVertexYBeforeCuts
x-position of the PV before the event selection.
TH1D * fHistoDCAZAfterCuts
xy-plane of the DCA after the track selection.
TH1D * fHistoDCAZBeforeCuts
xy-plane of the DCA before the track selection.
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.
TH1D * fHistoNumberOfTracksAfterEventCuts
Number of tracks in the event before any selection.
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)