AliPhysics  c2ade29 (c2ade29)
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 multiparticle correlations for the flow //
19 // harmonics v_1 to v_6. This version of the script compute the 2-, 4- and 6- particle //
20 // correlations for all the useful combinations of these six harmonics. It can take //
21 // Monte Carlo simulations data (e.g. HIJING), as well as the experimental Pb-Pb data //
22 // taken by the ALICE experiment. //
23 // The method used to compute the multiparticle correlations is the Generic Framework //
24 // based on Q-vectors. A setter lets open the possibility to cross-check the results //
25 // with nested loops. //
26 // //
27 // Author: Cindy Mordasini (cindy.mordasini@cern.ch) //
28 // Version: 27.02.2019 //
29 //--------------------------------------------------------------------------------------//
30 
32 #include "AliLog.h"
33 #include "AliAnalysisManager.h"
34 #include "Riostream.h"
35 #include "AliAODEvent.h"
36 #include "AliAODInputHandler.h"
37 #include "AliMCEvent.h"
38 #include "AliMCEventHandler.h"
39 #include "AliMultSelection.h"
40 #include "AliAODVertex.h"
41 #include "AliAODTrack.h"
42 #include "AliMCVertex.h"
43 #include "AliAODMCParticle.h"
44 #include "TMath.h"
45 #include "TComplex.h"
46 #include "TFile.h"
47 #include "TList.h"
48 #include <vector>
49 #include "TH1D.h"
50 #include "TH1I.h"
51 
52 using std::cout;
53 using std::endl;
54 
56 
57 //######################################################################################//
58 // Mandatory methods for AliAnalysisTaskSE.
59 //======================================================================================//
62 // Structure of the output file.
63  fMainList(NULL),
64  fQAListBeforeSelection(NULL),
65  fQAListAfterSelection(NULL),
66  fListCorrelations(NULL),
67 // General parameters.
68  fMaxNumberOfParticlesInCorrelations(8),
69  fHighestFlowHarmonic(6),
70  fUseParticleWeights(kFALSE),
71  fCrossCheckFourParticleCorrelations(kFALSE),
72  fCrossCheckWithNestedLoops(kFALSE),
73 // Type of files used in the analysis.
74  fProcessOnlyAOD(kFALSE),
75  fProcessOnlyMC(kFALSE),
76  fProcessBothMCandAOD(kFALSE),
77 // Determination of the centrality.
78  fCentralityFromVZero(kFALSE),
79  fCentralityFromSPD(kFALSE),
80  fCentralityMin(0.),
81  fCentralityMax(100.),
82 // Event selection.
83  fCutOnVertexX(kFALSE),
84  fVertexMinX(-44.),
85  fVertexMaxX(-44.),
86  fCutOnVertexY(kFALSE),
87  fVertexMinY(-44.),
88  fVertexMaxY(-44.),
89  fCutOnVertexZ(kFALSE),
90  fVertexMinZ(-10.),
91  fVertexMaxZ(10.),
92  fNumberOfTracksMin(6),
93 // Track selection.
94  fPtMin(0.2),
95  fPtMax(5.),
96  fEtaMin(-0.8),
97  fEtaMax(0.8),
98  fFilter(128),
99  fNumberOfTPCMin(70),
100  fChiSquarePInTPCMin(0.1),
101  fChiSquarePInTPCMax(4.),
102  fDCAxyMax(3.2),
103  fDCAzMax(2.4),
104 // TH1D with the observables for the event selection.
105  fHistoCentrality(NULL),
106  fHistoInitialNumberOfTracks(NULL),
107  fHistoNumberOfTracksBeforeTrackSelection(NULL),
108  fHistoFinalNumberOfTracks(NULL),
109  fHNOTNumberOfBins(30000),
110  fHNOTMax(30000.),
111  fHistoVertexXBeforeSelection(NULL),
112  fHistoVertexXAfterSelection(NULL),
113  fHistoVertexYBeforeSelection(NULL),
114  fHistoVertexYAfterSelection(NULL),
115  fHistoVertexZBeforeSelection(NULL),
116  fHistoVertexZAfterSelection(NULL),
117 // TH1D with the observables for the track selection.
118  fHistoPtBeforeSelection(NULL),
119  fHistoPtAfterSelection(NULL),
120  fHistoEtaBeforeSelection(NULL),
121  fHistoEtaAfterSelection(NULL),
122  fHistoPhiBeforeSelection(NULL),
123  fHistoPhiAfterSelection(NULL),
124  fHistoTPCClustersBeforeSelection(NULL),
125  fHistoTPCClustersAfterSelection(NULL),
126  fHistoTPCChiSquareBeforeSelection(NULL),
127  fHistoTPCChiSquareAfterSelection(NULL),
128  fHistoDCAXYBeforeSelection(NULL),
129  fHistoDCAXYAfterSelection(NULL),
130  fHistoDCAZBeforeSelection(NULL),
131  fHistoDCAZAfterSelection(NULL),
132 // TProfiles with the final multiparticle correlations.
133  fProfileTwoParticleCorrelations(NULL),
134  fProfileFourParticleCorrelations(NULL),
135  fProfileFourParticleCorrelationsCrossCheck(NULL),
136  fProfileSixParticleCorrelations(NULL),
137  fProfileTwoParticleCorrelationsNestedLoops(NULL),
138  fProfileFourParticleCorrelationsNestedLoops(NULL)
139 {
140 /* Dummy constructor of the class. */
141  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
142 
143 // Initialise 'fQvectors' to zero.
144  InitialiseArraysOfQvectors();
145 }
146 
147 //======================================================================================//
149  AliAnalysisTaskSE(name),
150 // Structure of the output file.
151  fMainList(NULL),
152  fQAListBeforeSelection(NULL),
153  fQAListAfterSelection(NULL),
154  fListCorrelations(NULL),
155 // General parameters.
156  fMaxNumberOfParticlesInCorrelations(8),
157  fHighestFlowHarmonic(6),
158  fUseParticleWeights(kFALSE),
159  fCrossCheckFourParticleCorrelations(kFALSE),
160  fCrossCheckWithNestedLoops(kFALSE),
161 // Type of files used in the analysis.
162  fProcessOnlyAOD(kFALSE),
163  fProcessOnlyMC(kFALSE),
164  fProcessBothMCandAOD(kFALSE),
165 // Determination of the centrality.
166  fCentralityFromVZero(kFALSE),
167  fCentralityFromSPD(kFALSE),
168  fCentralityMin(0.),
169  fCentralityMax(100.),
170 // Event selection.
171  fCutOnVertexX(kFALSE),
172  fVertexMinX(-44.),
173  fVertexMaxX(-44.),
174  fCutOnVertexY(kFALSE),
175  fVertexMinY(-44.),
176  fVertexMaxY(-44.),
177  fCutOnVertexZ(kFALSE),
178  fVertexMinZ(-10.),
179  fVertexMaxZ(10.),
180  fNumberOfTracksMin(6),
181 // Track selection.
182  fPtMin(0.2),
183  fPtMax(5.),
184  fEtaMin(-0.8),
185  fEtaMax(0.8),
186  fFilter(128),
187  fNumberOfTPCMin(70),
188  fChiSquarePInTPCMin(0.1),
189  fChiSquarePInTPCMax(4.),
190  fDCAxyMax(3.2),
191  fDCAzMax(2.4),
192 // TH1D with the observables for the event selection.
193  fHistoCentrality(NULL),
194  fHistoInitialNumberOfTracks(NULL),
195  fHistoNumberOfTracksBeforeTrackSelection(NULL),
196  fHistoFinalNumberOfTracks(NULL),
197  fHNOTNumberOfBins(30000),
198  fHNOTMax(30000.),
199  fHistoVertexXBeforeSelection(NULL),
200  fHistoVertexXAfterSelection(NULL),
201  fHistoVertexYBeforeSelection(NULL),
202  fHistoVertexYAfterSelection(NULL),
203  fHistoVertexZBeforeSelection(NULL),
204  fHistoVertexZAfterSelection(NULL),
205 // TH1D with the observables for the track selection.
206  fHistoPtBeforeSelection(NULL),
207  fHistoPtAfterSelection(NULL),
208  fHistoEtaBeforeSelection(NULL),
209  fHistoEtaAfterSelection(NULL),
210  fHistoPhiBeforeSelection(NULL),
211  fHistoPhiAfterSelection(NULL),
212  fHistoTPCClustersBeforeSelection(NULL),
213  fHistoTPCClustersAfterSelection(NULL),
214  fHistoTPCChiSquareBeforeSelection(NULL),
215  fHistoTPCChiSquareAfterSelection(NULL),
216  fHistoDCAXYBeforeSelection(NULL),
217  fHistoDCAXYAfterSelection(NULL),
218  fHistoDCAZBeforeSelection(NULL),
219  fHistoDCAZAfterSelection(NULL),
220 // TProfiles with the final multiparticle correlations.
221  fProfileTwoParticleCorrelations(NULL),
222  fProfileFourParticleCorrelations(NULL),
223  fProfileFourParticleCorrelationsCrossCheck(NULL),
224  fProfileSixParticleCorrelations(NULL),
225  fProfileTwoParticleCorrelationsNestedLoops(NULL),
226  fProfileFourParticleCorrelationsNestedLoops(NULL)
227 {
228 /* Constructor of the class. */
229  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
230 
231 // Create the mother list.
232  fMainList = new TList();
233  fMainList->SetName("outputAnalysis");
234  fMainList->SetOwner(kTRUE); // Gives ownership of the elements inside the TList to the TList itself.
235 
236 // Define the input and output slots.
237  DefineOutput(1, TList::Class());
238 
239 // Initialise 'fQvectors' to zero.
241 
242 // Define the procedure to follow if non-unit particle weights are used.
243  /*if (fUseParticleWeights)
244  {
245  // TBA, needed for data periods after 2010.
246  }*/
247 }
248 
249 //======================================================================================//
251 {
252 /* Destructor of the class. */
253  if (fMainList) {delete fMainList;}
254 }
255 
256 //======================================================================================//
258 {
259 /* Define the outputs of the task at the beginning of the analysis. */
260 // Avoid name clashes.
261  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
262  TH1::AddDirectory(kFALSE);
263 
264 // Book all the lists.
265  this->BookAllLists();
266 
267 // Book the histograms in all the daughter lists.
269  this->BookQAListAfterSelection();
270  this->BookListCorrelations();
271 
272 // Continue to avoid name clashes.
273  TH1::AddDirectory(oldHistAddStatus);
274  PostData(1, fMainList);
275 }
276 
277 //======================================================================================//
279 {
280 /* Execute the analysis for each event. */
282  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::UserExec(Option_t *)";
283 
284 // Select the type of file for the analysis (MC/AOD) from TaskSE.
285  AliAODEvent *currentAODEvent = dynamic_cast<AliAODEvent*>(InputEvent()); // Pointer to an AOD event.
286  AliMCEvent *currentMCEvent = MCEvent(); // Pointer to a Monte Carlo event.
287 
289  {
290  Fatal(sMethodName.Data(), "ERROR: only one fProcess must be kTRUE in 'SetAnalysisType'.");
291  }
292  else if (fProcessOnlyAOD) {AnalyseAODevent(currentAODEvent);}
293  else if (fProcessOnlyMC) {AnalyseMCevent(currentMCEvent);}
294  else if (fProcessBothMCandAOD) {Fatal(sMethodName.Data(),"ERROR: TBA.");}
295 
296 // PostData.
297  PostData(1, fMainList);
298 }
299 
300 //======================================================================================//
302 {
303 /* Save the outputs after the running over the events. */
304 // Access the mother list.
305  fMainList = (TList*)GetOutputData(1);
306  if (!fMainList) {exit(1);}
307 
308 // Create the output file and save the mother list inside.
309  TFile *outputFile = new TFile("AnalysisResults.root", "RECREATE");
310  fMainList->Write(fMainList->GetName(),TObject::kSingleKey);
311  delete outputFile;
312 }
313 
314 //######################################################################################//
315 // Methods called in the constructors.
316 //======================================================================================//
318 {
319 /* Initialise all the elements in 'fQvectors' to zero. */
320  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
321  {
322  for (Int_t iPower = 0; iPower < 9; iPower++)
323  {
324  fQvectors[iHarmo][iPower] = TComplex(0.,0.);
325  }
326  }
327 }
328 
329 //######################################################################################//
330 // Methods called in 'UserCreateOutputObjects'.
331 //======================================================================================//
333 {
334 /* Book all the lists in the output file */
335 // Check if the mother list exists.
336  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::BookAllLists()";
337  if (!fMainList) {Fatal(sMethodName.Data(), "Error: 'fMainList' is NULL.");}
338 
339 // Daughter list with the observables involved in the event selection.
341  fQAListBeforeSelection->SetName("fQAListBeforeSelection");
342  fQAListBeforeSelection->SetOwner(kTRUE);
343  fMainList->Add(fQAListBeforeSelection);
344 
345 // Daughter list with the observables involved in the track selection.
347  fQAListAfterSelection->SetName("fQAListAfterSelection");
348  fQAListAfterSelection->SetOwner(kTRUE);
350 
351 // Daughter list with the multiparticle correlations.
352  fListCorrelations = new TList();
353  fListCorrelations->SetName("fListCorrelations");
354  fListCorrelations->SetOwner(kTRUE);
356 }
357 
358 //======================================================================================//
360 {
361 /* Book the TH1D with the observables involved in the event selection. */
362 // Centrality distribution.
363  fHistoCentrality = new TH1D("fHistoCentrality", "Distribution of the centrality before the event selection", 100, 0., 100.);
364  fHistoCentrality->SetStats(kTRUE);
365  fHistoCentrality->GetXaxis()->SetTitle("Centrality percentile");
367 
368 // Distribution of the initial number of tracks.
369  fHistoInitialNumberOfTracks = new TH1I("fHistoInitialNumberOfTracks", "Initial distribution of the number of tracks", fHNOTNumberOfBins, 0., fHNOTMax);
370  fHistoInitialNumberOfTracks->SetStats(kTRUE);
371  fHistoInitialNumberOfTracks->GetXaxis()->SetTitle("Number of tracks");
373 
374 // Distribution of the number of tracks before the track selection.
375  fHistoNumberOfTracksBeforeTrackSelection = new TH1I("fHistoNumberOfTracksBeforeTrackSelection", "Distribution of the number of tracks before the track selection", fHNOTNumberOfBins, 0., fHNOTMax);
377  fHistoNumberOfTracksBeforeTrackSelection->GetXaxis()->SetTitle("Number of tracks");
379 
380 // x-position of the PV distribution.
381  fHistoVertexXBeforeSelection = new TH1D("fHistoVertexXBeforeSelection", "Distribution of PV_{x} before the selection", 1000, -20., 20.);
382  fHistoVertexXBeforeSelection->SetStats(kTRUE);
383  fHistoVertexXBeforeSelection->GetXaxis()->SetTitle("PV_{x}");
385 
386 // y-position of the PV distribution.
387  fHistoVertexYBeforeSelection = new TH1D("fHistoVertexYBeforeSelection", "Distribution of PV_{y} before the selection", 1000, -20., 20.);
388  fHistoVertexYBeforeSelection->SetStats(kTRUE);
389  fHistoVertexYBeforeSelection->GetXaxis()->SetTitle("PV_{y}");
391 
392 // z-position of the PV distribution.
393  fHistoVertexZBeforeSelection = new TH1D("fHistoVertexZBeforeSelection", "Distribution of PV_{z} before the selection", 1000, -20., 20.);
394  fHistoVertexZBeforeSelection->SetStats(kTRUE);
395  fHistoVertexZBeforeSelection->GetXaxis()->SetTitle("PV_{z}");
397 
398 // Transverse momentum distribution.
399  fHistoPtBeforeSelection = new TH1D("fHistoPtBeforeSelection", "Distribution of p_{T} before the track selection", 1000, 0., 20.);
400  fHistoPtBeforeSelection->SetStats(kTRUE);
401  fHistoPtBeforeSelection->GetXaxis()->SetTitle("p_{T}");
403 
404 // Pseudorapidity distribution.
405  fHistoEtaBeforeSelection = new TH1D("fHistoEtaBeforeSelection", "Distribution of #eta before the track selection", 1000, -1., 1.);
406  fHistoEtaBeforeSelection->SetStats(kTRUE);
407  fHistoEtaBeforeSelection->GetXaxis()->SetTitle("#eta");
409 
410 // Azimuthal angles distribution.
411  fHistoPhiBeforeSelection = new TH1D("fHistoPhiBeforeSelection", "Distributiion of #phi before the track selection", 1000, 0., 6.3);
412  fHistoPhiBeforeSelection->SetStats(kTRUE);
413  fHistoPhiBeforeSelection->GetXaxis()->SetTitle("#phi");
415 
416 // Distribution of the number of TPC clusters.
417  fHistoTPCClustersBeforeSelection = new TH1I("fHistoTPCClustersBeforeSelection", "Distribution of the number of TPC clusters before the track selection", 1000, 0., 170.);
418  fHistoTPCClustersBeforeSelection->SetStats(kTRUE);
419  fHistoTPCClustersBeforeSelection->GetXaxis()->SetTitle("Number of TPC clusters");
421 
422 // Distribution of the chi square of the track momentum in the TPC.
423  fHistoTPCChiSquareBeforeSelection = new TH1D("fHistoTPCChiSquareBeforeSelection", "Distribution of the #chi^{2} of the track momentum in the TPC before the track selection", 1000, 0., 20.);
424  fHistoTPCChiSquareBeforeSelection->SetStats(kTRUE);
425  fHistoTPCChiSquareBeforeSelection->GetXaxis()->SetTitle("#chi^{2}/NDF in TPC");
427 
428 // xy-plane of the DCA distribution.
429  fHistoDCAXYBeforeSelection = new TH1D("fHistoDCAXYBeforeSelection", "Distribution of DCA_{xy} before the track selection", 1000, 0., 10.);
430  fHistoDCAXYBeforeSelection->SetStats(kTRUE);
431  fHistoDCAXYBeforeSelection->GetXaxis()->SetTitle("DCA_{xy}");
433 
434 // z-coordinate of the DCA distribution.
435  fHistoDCAZBeforeSelection = new TH1D("fHistoDCAZBeforeSelection", "Distribution of DCA_{z} before the track selection", 1000, 0., 10.);
436  fHistoDCAZBeforeSelection->SetStats(kTRUE);
437  fHistoDCAZBeforeSelection->GetXaxis()->SetTitle("DCA_{z}");
439 }
440 
441 //======================================================================================//
443 {
444 /* Book the TH1D with the observables involved in the track selection. */
445 // Final number of tracks after both the event and track selection.
446  fHistoFinalNumberOfTracks = new TH1I("fHistoFinalNumberOfTracks", "Final distribution of the number of tracks", fHNOTNumberOfBins, 0., fHNOTMax);
447  fHistoFinalNumberOfTracks->SetStats(kTRUE);
448  fHistoFinalNumberOfTracks->GetXaxis()->SetTitle("Number of tracks");
450 
451 // x-position of the PV distribution.
452  fHistoVertexXAfterSelection = new TH1D("fHistoVertexXAfterSelection", "Distribution of PV_{x} after the full selection", 1000, -20., 20.);
453  fHistoVertexXAfterSelection->SetStats(kTRUE);
454  fHistoVertexXAfterSelection->GetXaxis()->SetTitle("PV_{x}");
456 
457 // y-position of the PV distribution.
458  fHistoVertexYAfterSelection = new TH1D("fHistoVertexYAfterSelection", "Distribution of PV_{y} after the full selection", 1000, -20., 20.);
459  fHistoVertexYAfterSelection->SetStats(kTRUE);
460  fHistoVertexYAfterSelection->GetXaxis()->SetTitle("PV_{y}");
462 
463 // z-position of the PV distribution.
464  fHistoVertexZAfterSelection = new TH1D("fHistoVertexZAfterSelection", "Distribution of PV_{z} after the full selection", 1000, -20., 20.);
465  fHistoVertexZAfterSelection->SetStats(kTRUE);
466  fHistoVertexZAfterSelection->GetXaxis()->SetTitle("PV_{z}");
468 
469 // Transverse momentum distribution.
470  fHistoPtAfterSelection = new TH1D("fHistoPtAfterSelection", "Distribution of p_{T} after the full selection", 1000, 0., 20.);
471  fHistoPtAfterSelection->SetStats(kTRUE);
472  fHistoPtAfterSelection->GetXaxis()->SetTitle("p_{T}");
474 
475 // Pseudorapidity distribution.
476  fHistoEtaAfterSelection = new TH1D("fHistoEtaAfterSelection", "Distribution of #eta after the full selection", 1000, -1., 1.);
477  fHistoEtaAfterSelection->SetStats(kTRUE);
478  fHistoEtaAfterSelection->GetXaxis()->SetTitle("#eta");
480 
481 // Azimuthal angles distribution.
482  fHistoPhiAfterSelection = new TH1D("fHistoPhiAfterSelection", "Distribution of #phi after the full selection", 1000, 0., 6.3);
483  fHistoPhiAfterSelection->SetStats(kTRUE);
484  fHistoPhiAfterSelection->GetXaxis()->SetTitle("#phi");
486 
487 // Distribution of the number of TPC clusters.
488  fHistoTPCClustersAfterSelection = new TH1I("fHistoTPCClustersAfterSelection", "Distribution of the number of TPC clusters after the full selection", 1000, 0., 170.);
489  fHistoTPCClustersAfterSelection->SetStats(kTRUE);
490  fHistoTPCClustersAfterSelection->GetXaxis()->SetTitle("Number of TPC clusters");
492 
493 // Distribution of the chi square of the track momentum in the TPC.
494  fHistoTPCChiSquareAfterSelection = new TH1D("fHistoTPCChiSquareAfterSelection", "Distribution of the #chi^{2} of the track momentum in the TPC after the full selection", 1000, 0., 20.);
495  fHistoTPCChiSquareAfterSelection->SetStats(kTRUE);
496  fHistoTPCChiSquareAfterSelection->GetXaxis()->SetTitle("#chi^{2}/NDF in TPC");
498 
499 // xy-plane of the DCA distribution.
500  fHistoDCAXYAfterSelection = new TH1D("fHistoDCAXYAfterSelection", "Distribution of DCA_{xy} after the full selection", 1000, 0., 10.);
501  fHistoDCAXYAfterSelection->SetStats(kTRUE);
502  fHistoDCAXYAfterSelection->GetXaxis()->SetTitle("DCA_{xy}");
504 
505 // z-coordinate of the DCA distribution.
506  fHistoDCAZAfterSelection = new TH1D("fHistoDCAZAfterSelection", "Distribution of DCA_{z} after the full selection", 1000, 0., 10.);
507  fHistoDCAZAfterSelection->SetStats(kTRUE);
508  fHistoDCAZAfterSelection->GetXaxis()->SetTitle("DCA_{z}");
510 }
511 
512 //======================================================================================//
514 {
515 /* Book the TProfiles with the multiparticle correlations. */
516 // 2-particle correlations.
517  fProfileTwoParticleCorrelations = new TProfile("fProfileTwoParticleCorrelations", "2-particle correlations", 6, 0., 6.);
518  fProfileTwoParticleCorrelations->SetStats(kTRUE);
520  fProfileTwoParticleCorrelations->GetXaxis()->SetTitle("n");
521  fProfileTwoParticleCorrelations->GetYaxis()->SetTitle("#LT#LT2#GT#GT_{n,-n}");
523 
525  {
526  fProfileTwoParticleCorrelationsNestedLoops = new TProfile("fProfileTwoParticleCorrelationsNestedLoops", "2-particle correlations with nested loops", 6, 0., 6.);
529  fProfileTwoParticleCorrelationsNestedLoops->GetXaxis()->SetTitle("n");
530  fProfileTwoParticleCorrelationsNestedLoops->GetYaxis()->SetTitle("#LT#LT2#GT#GT_{n,-n}");
532  }
533 
534 // 4-particle correlations.
535  fProfileFourParticleCorrelations = new TProfile("fProfileFourParticleCorrelations", "4-particle correlations", 21, 0., 21.);
536  fProfileFourParticleCorrelations->SetStats(kTRUE);
538  fProfileFourParticleCorrelations->GetXaxis()->SetTitle("(m,n)");
539  fProfileFourParticleCorrelations->GetYaxis()->SetTitle("#LT#LT4#GT#GT_{m,n,-m,-n}");
541 
543  {
544  fProfileFourParticleCorrelationsCrossCheck = new TProfile("fProfileFourParticleCorrelationsCrossCheck", "4-particle correlations for cross-check", 15, 0., 15.);
547  fProfileFourParticleCorrelationsCrossCheck->GetXaxis()->SetTitle("(m,n)");
548  fProfileFourParticleCorrelationsCrossCheck->GetYaxis()->SetTitle("#LT#LT4#GT#GT_{m,n,-m,-n}");
550  }
551 
553  {
554  fProfileFourParticleCorrelationsNestedLoops = new TProfile("fProfileFourParticleCorrelationsNestedLoops", "4-particle correlations with nested loops", 21, 0., 21.);
557  fProfileFourParticleCorrelationsNestedLoops->GetXaxis()->SetTitle("(m,n)");
558  fProfileFourParticleCorrelationsNestedLoops->GetYaxis()->SetTitle("#LT#LT4#GT#GT_{m,n,-m,-n}");
560  }
561 
562 // 6-particle correlations.
563  fProfileSixParticleCorrelations = new TProfile("fProfileSixParticleCorrelations", "6-particle correlations", 20, 0., 20.);
564  fProfileSixParticleCorrelations->SetStats(kTRUE);
566  fProfileSixParticleCorrelations->GetXaxis()->SetTitle("(l,m,n)");
567  fProfileSixParticleCorrelations->GetYaxis()->SetTitle("#LT#LT6#GT#GT_{l,m,n,-l,-m,-n}");
569 }
570 
571 //######################################################################################//
572 // Methods called in 'UserExec'.
573 //======================================================================================//
575 {
576 /* Execute the analysis for the provided AOD event. */
577  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::AnalyseAODevent(AliAODEvent *aAODevent)";
578 
579 // Check if there is an event or not.
580  if (!aAODevent) {Fatal(sMethodName.Data(), "ERROR: no AOD event found.");}
581 
582 // Select the detector to use for the estimation of the centrality.
583  TString centralityEstimator = "centralityEstimator"; // Name of the detector used for the centrality estimation.
585  {
586  Fatal(sMethodName.Data(), "ERROR: only one detector must be selected in 'SetCentralityEstimation'.");
587  }
588  else if (fCentralityFromVZero) {centralityEstimator = "V0M";}
589  else if (fCentralityFromSPD) {centralityEstimator = "CL1";}
590 
591 // Determine if the event belongs to this centrality range.
592  AliMultSelection *ams = (AliMultSelection*)aAODevent->FindListObject("MultSelection");
593  if (!ams) {return;} // Protection against NULL pointer.
594  Double_t aCentrality = ams->GetMultiplicityPercentile(Form("%s", centralityEstimator.Data())); // Centrality of the given event.
595  if ((aCentrality >= fCentralityMin) && (aCentrality < fCentralityMax))
596  {
597  fHistoCentrality->Fill(aCentrality);
598  }
599  else {return;} // This event does not belong to this centrality range.
600 
601 // Get the number of tracks before the event selection.
602  long long initialNumberOfTracks = aAODevent->GetNumberOfTracks();
603  fHistoInitialNumberOfTracks->Fill(initialNumberOfTracks);
604 
605 // Cuts on the position of the Primary Vertex.
606  AliAODVertex *avtx = (AliAODVertex*)aAODevent->GetPrimaryVertex(); // 3d position of the PV.
607  fHistoVertexXBeforeSelection->Fill(avtx->GetX());
608  fHistoVertexYBeforeSelection->Fill(avtx->GetY());
609  fHistoVertexZBeforeSelection->Fill(avtx->GetZ());
610 
611  if (fCutOnVertexX)
612  {
613  if ((avtx->GetX() < fVertexMinX) || (avtx->GetX() > fVertexMaxX)) {return;}
614  }
615  if (fCutOnVertexY)
616  {
617  if ((avtx->GetY() < fVertexMinY) || (avtx->GetY() > fVertexMaxY)) {return;}
618  }
619  if (fCutOnVertexZ)
620  {
621  if ((avtx->GetZ() < fVertexMinZ) || (avtx->GetZ() > fVertexMaxZ)) {return;}
622  }
623 
625 
626 // Preparations for the track selection.
627  long long numberOfTracksBeforeTrackSelection = aAODevent->GetNumberOfTracks(); // Number of tracks before the track selection.
628  fHistoNumberOfTracksBeforeTrackSelection->Fill(numberOfTracksBeforeTrackSelection);
629  long long finalNumberOfTracks = 0; // Number of tracks after the full selection.
630  Int_t *IsTrackSelected = new Int_t[numberOfTracksBeforeTrackSelection](); // Flag to indicate a track passed the track selection (1) or not (0).
631 
632  Double_t pT = 0.; // Transverse momentum.
633  Double_t eta = 0.; // Pseudorapidity.
634  Double_t phi = 0.; // Azimuthal angle.
635  Int_t numberOfTPCClusters = 0; // Number of TPC clusters.
636  Double_t chiSquareInTPC = 0.; // Chi square of the track momentum in the TPC.
637  Double_t DCAx = 0.; // x-value of the DCA.
638  Double_t DCAy = 0.; // y-value of the DCA.
639  Double_t DCAz = 0.; // z-value of the DCA.
640  Double_t DCAxy = 0.; // xy-value of the DCA.
641 
642 // Look at each track in the event to mark them as selected or not.
643  for (long long iTrack = 0; iTrack < numberOfTracksBeforeTrackSelection; iTrack++)
644  {
645  AliAODTrack *currentTrack = dynamic_cast<AliAODTrack*>(aAODevent->GetTrack(iTrack)); // Pointer to the AOD track.
646  if (!currentTrack) {continue;} // Protection against NULL pointer.
647  if (!currentTrack->TestFilterBit(fFilter)) {continue;} // Filter bit 128 denotes TPC-only tracks.
648 
649  // Get all the observables for the track selection.
650  pT = currentTrack->Pt();
651  eta = currentTrack->Eta();
652  phi = currentTrack->Phi();
653  numberOfTPCClusters = currentTrack->GetTPCNcls();
654  chiSquareInTPC = currentTrack->Chi2perNDF();
655  DCAx = currentTrack->XAtDCA();
656  DCAy = currentTrack->YAtDCA();
657  DCAz = currentTrack->ZAtDCA();
658 
659  DCAxy = TMath::Sqrt((DCAx*DCAx) + (DCAy*DCAy));
660 
661  // Fill the histograms before the track selection.
662  fHistoPtBeforeSelection->Fill(pT);
663  fHistoEtaBeforeSelection->Fill(eta);
664  fHistoPhiBeforeSelection->Fill(phi);
665  fHistoTPCClustersBeforeSelection->Fill(numberOfTPCClusters);
666  fHistoTPCChiSquareBeforeSelection->Fill(chiSquareInTPC);
667  fHistoDCAXYBeforeSelection->Fill(DCAxy);
668  fHistoDCAZBeforeSelection->Fill(DCAz);
669 
670  // Apply the track selection to the provided track.
671  if (ApplyTrackSelection(pT, eta, numberOfTPCClusters, chiSquareInTPC, DCAxy, DCAz)) // The track passed the selection.
672  {
673  IsTrackSelected[iTrack] = 1;
674  finalNumberOfTracks++;
675  }
676  else {IsTrackSelected[iTrack] = 0;} // The track failed the selection.
677  }
678 
679 // Remove the events with not enough tracks for the event weight.
680  if (finalNumberOfTracks <= fNumberOfTracksMin) {return;}
681 
682 // Fill all the event histograms after the full selection.
683  fHistoFinalNumberOfTracks->Fill(finalNumberOfTracks);
684  fHistoVertexXAfterSelection->Fill(avtx->GetX());
685  fHistoVertexYAfterSelection->Fill(avtx->GetY());
686  fHistoVertexZAfterSelection->Fill(avtx->GetZ());
687 
688 // Define the arrays for the azimuthal angles and particle weights to use in the analysis.
689  Double_t *phiArray = new Double_t[finalNumberOfTracks](); // Azimuthal angles.
690  Double_t *particleWeightArray = new Double_t[finalNumberOfTracks](); // Particle weights.
691  Int_t indexInNewArrays = 0; // New index of the track if it passed the selection.
692 
693 // Loop over the tracks to keep only the selected ones.
694  for (long long iTrack = 0; iTrack < numberOfTracksBeforeTrackSelection; iTrack++)
695  {
696  AliAODTrack *aTrack = dynamic_cast<AliAODTrack*>(aAODevent->GetTrack(iTrack)); // Pointer to the AOD track.
697  if (!aTrack) {continue;} // Protection against NULL pointer.
698  if (!aTrack->TestFilterBit(fFilter)) {continue;} // Filter bit 128 denotes TPC-only tracks.
699 
700  if (IsTrackSelected[iTrack] == 1) // The particle passed the selection.
701  {
702  // Get all the observables used in the track selection.
703  pT = aTrack->Pt();
704  eta = aTrack->Eta();
705  phiArray[indexInNewArrays] = aTrack->Phi();
706  numberOfTPCClusters = aTrack->GetTPCNcls();
707  chiSquareInTPC = aTrack->Chi2perNDF();
708  DCAx = aTrack->XAtDCA();
709  DCAy = aTrack->YAtDCA();
710  DCAz = aTrack->ZAtDCA();
711  DCAxy = TMath::Sqrt((DCAx*DCAx) + (DCAy*DCAy));
712 
713  if (fUseParticleWeights) {Fatal(sMethodName.Data(), "ERROR: TBA.");}
714  else {particleWeightArray[indexInNewArrays] = 1.;}
715 
716  // Fill all the track histograms after the full selection.
717  fHistoPtAfterSelection->Fill(pT);
718  fHistoEtaAfterSelection->Fill(eta);
719  fHistoPhiAfterSelection->Fill(phiArray[indexInNewArrays]);
720  fHistoTPCClustersAfterSelection->Fill(numberOfTPCClusters);
721  fHistoTPCChiSquareAfterSelection->Fill(chiSquareInTPC);
722  fHistoDCAXYAfterSelection->Fill(DCAxy);
723  fHistoDCAZAfterSelection->Fill(DCAz);
724 
725  // Increase the value of 'indexInNewArrays' by one.
726  indexInNewArrays++;
727  }
728  else {continue;}
729  }
730 
731 // Calculate the Q-vectors for the current event.
732  CalculateQvectors(finalNumberOfTracks, phiArray, particleWeightArray);
733 
734 // Compute all the multiparticle correlations for the current event.
735  ComputeMultiparticleCorrelations(finalNumberOfTracks, phiArray, particleWeightArray);
736 
737 // Reset everything to zero for the next event.
738  numberOfTracksBeforeTrackSelection = 0;
739  finalNumberOfTracks = 0;
740  delete [] IsTrackSelected;
741  pT = 0.;
742  eta = 0.;
743  phi = 0.;
744  numberOfTPCClusters = 0;
745  chiSquareInTPC = 0.;
746  DCAx = 0.;
747  DCAy = 0.;
748  DCAz = 0.;
749  delete [] phiArray;
750  delete [] particleWeightArray;
751  indexInNewArrays = 0;
752 }
753 
754 //======================================================================================//
756 {
757 /* Execute the analysis for the provided MC event. */
758  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::AnalyseMCevent(AliMCEvent *aMCevent)";
759 
760 // Check if there is an event or not.
761  if (!aMCevent) {Fatal(sMethodName.Data(), "ERROR: no MC event found.");}
762 
763 // Select the detector to use for the estimation of the centrality.
764  TString centralityEstimator = "centralityEstimator"; // Name of the detector used for the centrality estimation.
766  {
767  Fatal(sMethodName.Data(), "ERROR: only one detector must be selected in 'SetCentralityEstimation'.");
768  }
769  else if (fCentralityFromVZero) {centralityEstimator = "V0M";}
770  else if (fCentralityFromSPD) {centralityEstimator = "CL1";}
771 
772 // Determine if the event belongs to this centrality range (for reconstructed particles only).
773  if (!fProcessOnlyMC)
774  {
775  AliMultSelection *ams = (AliMultSelection*)aMCevent->FindListObject("MultSelection");
776  if (!ams) {return;} // Protection against NULL pointer.
777  Double_t aCentrality = ams->GetMultiplicityPercentile(Form("%s", centralityEstimator.Data())); // Centrality of the given event.
778  if ((aCentrality >= fCentralityMin) && (aCentrality < fCentralityMax))
779  {
780  fHistoCentrality->Fill(aCentrality);
781  }
782  else {return;} // This event does not belong to this centrality range.
783  }
784 
785 // Get the number of tracks before the event selection.
786  long long initialNumberOfTracks = aMCevent->GetNumberOfTracks();
787  fHistoInitialNumberOfTracks->Fill(initialNumberOfTracks);
788 
789 // Cuts on the position of the Primary Vertex.
790  AliMCVertex *avtx = (AliMCVertex*)aMCevent->GetPrimaryVertex(); // 3d position of the PV.
791  fHistoVertexXBeforeSelection->Fill(avtx->GetX());
792  fHistoVertexYBeforeSelection->Fill(avtx->GetY());
793  fHistoVertexZBeforeSelection->Fill(avtx->GetZ());
794 
795  if (fCutOnVertexX)
796  {
797  if ((avtx->GetX() < fVertexMinX) || (avtx->GetX() > fVertexMaxX)) {return;}
798  }
799  if (fCutOnVertexY)
800  {
801  if ((avtx->GetY() < fVertexMinY) || (avtx->GetY() > fVertexMaxY)) {return;}
802  }
803  if (fCutOnVertexZ)
804  {
805  if ((avtx->GetZ() < fVertexMinZ) || (avtx->GetZ() > fVertexMaxZ)) {return;}
806  }
807 
809 
810 // Preparations for the track selection.
811  long long numberOfTracksBeforeTrackSelection = aMCevent->GetNumberOfTracks(); // Number of tracks before the track selection.
812  fHistoNumberOfTracksBeforeTrackSelection->Fill(numberOfTracksBeforeTrackSelection);
813  long long finalNumberOfTracks = 0; // Number of tracks after the full selection.
814  Int_t *IsTrackSelected = new Int_t[numberOfTracksBeforeTrackSelection](); // Flag to indicate a track passed the track selection (1) or not (0).
815 
816  Double_t pT = 0.; // Transverse momentum.
817  Double_t eta = 0.; // Pseudorapidity.
818  Double_t phi = 0.; // Azimuthal angle.
819 
820 // Look at each track in the event to mark them as selected or not.
821  for (long long iTrack = 0; iTrack < numberOfTracksBeforeTrackSelection; iTrack++)
822  {
823  AliAODMCParticle *currentTrack = dynamic_cast<AliAODMCParticle*>(aMCevent->GetTrack(iTrack)); // Pointer to the MC track.
824  if (!currentTrack) {continue;} // Protection against NULL pointer.
825 
826  // Get all the observables for the track selection.
827  pT = currentTrack->Pt();
828  eta = currentTrack->Eta();
829  phi = currentTrack->Phi();
830 
831  // Fill the histograms before the track selection.
832  fHistoPtBeforeSelection->Fill(pT);
833  fHistoEtaBeforeSelection->Fill(eta);
834  fHistoPhiBeforeSelection->Fill(phi);
835 
836 
837  // Apply the track selection to the provided track.
838  if ((fPtMin <= pT) && (pT <= fPtMax) && (fEtaMin <= eta) && (eta <= fEtaMax)) // The track passed the selection.
839  {
840  IsTrackSelected[iTrack] = 1;
841  finalNumberOfTracks++;
842  }
843  else {IsTrackSelected[iTrack] = 0;} // The track failed the selection.
844  }
845 
846 // Remove the events with not enough tracks for the event weight.
847  if (finalNumberOfTracks <= fNumberOfTracksMin) {return;}
848 
849 // Fill all the event histograms after the full selection.
850  fHistoFinalNumberOfTracks->Fill(finalNumberOfTracks);
851  fHistoVertexXAfterSelection->Fill(avtx->GetX());
852  fHistoVertexYAfterSelection->Fill(avtx->GetY());
853  fHistoVertexZAfterSelection->Fill(avtx->GetZ());
854 
855 // Define the arrays for the azimuthal angles and particle weights to use in the analysis.
856  Double_t *phiArray = new Double_t[finalNumberOfTracks](); // Azimuthal angles.
857  Double_t *particleWeightArray = new Double_t[finalNumberOfTracks](); // Particle weights.
858  Int_t indexInNewArrays = 0; // New index of the track if it passed the selection.
859 
860 // Loop over the tracks to keep only the selected ones.
861  for (long long iTrack = 0; iTrack < numberOfTracksBeforeTrackSelection; iTrack++)
862  {
863  AliAODMCParticle *aTrack = dynamic_cast<AliAODMCParticle*>(aMCevent->GetTrack(iTrack)); // Pointer to the MC track.
864  if (!aTrack) {continue;} // Protection against NULL pointer.
865 
866  if (IsTrackSelected[iTrack] == 1) // The particle passed the selection.
867  {
868  // Get all the observables used in the track selection.
869  pT = aTrack->Pt();
870  eta = aTrack->Eta();
871  phiArray[indexInNewArrays] = aTrack->Phi();
872 
873  if (fUseParticleWeights) {Fatal(sMethodName.Data(), "ERROR: TBA.");}
874  else {particleWeightArray[indexInNewArrays] = 1.;}
875 
876  // Fill all the track histograms after the full selection.
877  fHistoPtAfterSelection->Fill(pT);
878  fHistoEtaAfterSelection->Fill(eta);
879  fHistoPhiAfterSelection->Fill(phiArray[indexInNewArrays]);
880 
881  // Increase the value of 'indexInNewArrays' by one.
882  indexInNewArrays++;
883  }
884  else {continue;}
885  }
886 
887 // Calculate the Q-vectors for the current event.
888  CalculateQvectors(finalNumberOfTracks, phiArray, particleWeightArray);
889 
890 // Compute all the multiparticle correlations for the current event.
891  ComputeMultiparticleCorrelations(finalNumberOfTracks, phiArray, particleWeightArray);
892 
893 // Reset everything to zero for the next event.
894  numberOfTracksBeforeTrackSelection = 0;
895  finalNumberOfTracks = 0;
896  delete [] IsTrackSelected;
897  pT = 0.;
898  eta = 0.;
899  phi = 0.;
900  delete [] phiArray;
901  delete [] particleWeightArray;
902  indexInNewArrays = 0;
903 }
904 
905 //======================================================================================//
906 Bool_t AliAnalysisTaskTwoMultiCorrelations::ApplyTrackSelection(Double_t momentum, Double_t pseudorapidity, Int_t NclustersInTPC, Double_t TPCchiSquare, Double_t xyDCA, Double_t zDCA)
907 {
908 /* Apply the track selection to the arguments and return if it is passed or not. */
909  Bool_t cutOnPt = (fPtMin <= momentum) && (momentum <= fPtMax);
910  Bool_t cutOnEta = (fEtaMin <= pseudorapidity) && (pseudorapidity <= fEtaMax);
911  Bool_t cutOnNumberOfTPC = (fNumberOfTPCMin < NclustersInTPC);
912  Bool_t cutOnChiSquare = (fChiSquarePInTPCMin <= TPCchiSquare) && (TPCchiSquare <= fChiSquarePInTPCMax);
913  Bool_t cutOnDCAxy = (xyDCA < fDCAxyMax);
914  Bool_t cutOnDCAz = (zDCA < fDCAzMax);
915 
916  return cutOnPt && cutOnEta && cutOnNumberOfTPC && cutOnChiSquare && cutOnDCAxy && cutOnDCAz;
917 }
918 
919 //======================================================================================//
920 void AliAnalysisTaskTwoMultiCorrelations::CalculateQvectors(long long numberOfParticles, Double_t angles[], Double_t pWeights[])
921 {
922 /* Calculate all the Q-vectors for the given arrays of azimuthal angles and particle weights. */
923  Double_t iAngle = 0.; // Azimuthal angle of the current particle.
924  Double_t iWeight = 0.; // Particle weight of the current particle.
925  Double_t iWeightToPowerP = 0.; // Particle weight rised to the power p.
926 
927 // Ensure all the Q-vectors are initially zero.
928  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
929  {
930  for (Int_t iPower = 0; iPower < 9; iPower++)
931  {
932  fQvectors[iHarmo][iPower] = TComplex(0.,0.);
933  }
934  }
935 
936 // Compute the Q-vectors.
937  for (long long iTrack = 0; iTrack < numberOfParticles; iTrack++)
938  {
939  iAngle = angles[iTrack];
940  iWeight = pWeights[iTrack];
941  for (Int_t iHarmo = 0; iHarmo < 49; iHarmo++)
942  {
943  for (Int_t iPower = 0; iPower < 9; iPower++)
944  {
945  iWeightToPowerP = TMath::Power(iWeight, iPower);
946  fQvectors[iHarmo][iPower] += TComplex(iWeightToPowerP*TMath::Cos(iHarmo*iAngle), iWeightToPowerP*TMath::Sin(iHarmo*iAngle));
947  }
948  }
949  }
950 }
951 
952 //======================================================================================//
954 {
955 /* Simplify the use of the Q-vectors. */
956  if (n >= 0) {return fQvectors[n][p];}
957  return TComplex::Conjugate(fQvectors[-n][p]); // Use that Q*(n,p) = Q(-n,p).
958 }
959 
960 //======================================================================================//
961 void AliAnalysisTaskTwoMultiCorrelations::ComputeMultiparticleCorrelations(long long numberOfParticles, Double_t angles[], Double_t pWeights[])
962 {
963 /* Compute the considered 2-, 4- and 6-particle correlations for the harmonics v_1 to v_6. */
964 // Compute the 2-particle correlations.
965  Int_t twoZerosArray[2] = {0}; // Zero array for the denominator.
966  Double_t twoParticleDenominator = CalculateRecursion(2, twoZerosArray).Re(); // 2-particle event weight.
967 
968  Int_t twoParticleArray[2] = {0}; // Array to save the harmonics with the format: cos(nphi1 - nphi2).
969  TComplex twoParticleComplex = TComplex(0., 0.); // Complex value for the 2-particle correlation.
970  Double_t twoParticleValue = 0.; // Real part of the 2-particle correlation.
971  Double_t iBinTwoMiddle = 0.; // Index of the corresponding bin in the TProfile.
972  for (Int_t n = 1; n < 7; n++)
973  {
974  // Compute the correlation.
975  twoParticleArray[0] = n;
976  twoParticleArray[1] = -1*n;
977  twoParticleComplex = (CalculateRecursion(2, twoParticleArray))/twoParticleDenominator;
978  twoParticleValue = twoParticleComplex.Re();
979 
980  // Fill the corresponding bin in the TProfile.
981  iBinTwoMiddle = (1.*n) - 0.5;
982  fProfileTwoParticleCorrelations->Fill(iBinTwoMiddle, twoParticleValue, twoParticleDenominator);
983  fProfileTwoParticleCorrelations->GetXaxis()->SetBinLabel(n, Form("%d", n));
984 
985  // Cross-check with nested loops if needed.
987  {
988  // Fill the corresponding bin in the TProfile.
989  ComputeTwoNestedLoops(numberOfParticles, twoParticleArray, angles, pWeights, fProfileTwoParticleCorrelationsNestedLoops, iBinTwoMiddle);
990  fProfileTwoParticleCorrelationsNestedLoops->GetXaxis()->SetBinLabel(n, Form("%d", n));
991  }
992 
993  // Reset of the variables for the next harmonic.
994  twoParticleComplex = TComplex(0., 0.);
995  twoParticleValue = 0.;
996  iBinTwoMiddle = 0.;
997  }
998 
999 // Compute the 4-particle correlations.
1000  Int_t fourZerosArray[4] = {0}; // Zero array for the denominator.
1001  Double_t fourParticleDenominator = CalculateRecursion(4, fourZerosArray).Re(); // 4-particle event weight.
1002 
1003  Int_t fourParticleArray[4] = {0}; // Array to save the harmonics with the format: cos(mphi1 + nphi2 - mphi3 - nphi4).
1004  TComplex fourParticleComplex = TComplex(0., 0.); // Complex value for the 4-particle correlation.
1005  Double_t fourParticleValue = 0.; // Real part of the 4-particle correlation.
1006  Int_t iBinFour = 1; // Index of the corresponding bin in the TProfile.
1007  Int_t iBinFourMiddle = 0.; // Middle of the corresponding bin in the TProfile.
1008  for(Int_t m = 1; m < 7; m++)
1009  {
1010  for (Int_t n = m; n < 7; n++)
1011  {
1012  // Compute the correlation.
1013  fourParticleArray[0] = m;
1014  fourParticleArray[1] = n;
1015  fourParticleArray[2] = -1*m;
1016  fourParticleArray[3] = -1*n;
1017  fourParticleComplex = (CalculateRecursion(4, fourParticleArray))/fourParticleDenominator;
1018  fourParticleValue = fourParticleComplex.Re();
1019 
1020  // Fill the corresponding bin in the TProfile.
1021  iBinFourMiddle = (1.*iBinFour) - 0.5;
1022  fProfileFourParticleCorrelations->Fill(iBinFourMiddle, fourParticleValue, fourParticleDenominator);
1023  fProfileFourParticleCorrelations->GetXaxis()->SetBinLabel(iBinFour, Form("(%d,%d)", m, n));
1024 
1025  // Cross-check with nested loops if needed.
1027  {
1028  // Fill the corresponding bin in the TProfile.
1029  ComputeFourNestedLoops(numberOfParticles, fourParticleArray, angles, pWeights, fProfileFourParticleCorrelationsNestedLoops, iBinFourMiddle);
1030  fProfileFourParticleCorrelationsNestedLoops->GetXaxis()->SetBinLabel(iBinFour, Form("(%d,%d)", m, n));
1031  }
1032 
1033  // Reset of the variables for the next pair of harmonics.
1034  iBinFour++;
1035  fourParticleComplex = TComplex(0., 0.);
1036  fourParticleValue = 0.;
1037  iBinFourMiddle = 0.;
1038  }
1039  }
1040 
1041 // Compute the 4-particle correlations for the cross-check if needed.
1043  {
1044  Int_t fourParticleCrossCheckArray[4] = {0}; // Array to save the harmonics with the format: cos(mphi1 + nphi2 - mphi3 - nphi4).
1045  TComplex fourParticleCrossCheckComplex = TComplex(0., 0.); // Complex value for the 4-particle correlation.
1046  Double_t fourParticleCrossCheckValue = 0.; // Real part of the 4-particle correlation.
1047  Int_t iBinFourCrossCheck = 1; // Index of the corresponding bin in the TProfile.
1048  Int_t iBinFourCrossCheckMiddle = 0.; // Middle of the corresponding bin in the TProfile.
1049  for(Int_t m = 2; m < 7; m++)
1050  {
1051  for (Int_t n = 1; n < m; n++)
1052  {
1053  // Compute the correlation.
1054  fourParticleCrossCheckArray[0] = m;
1055  fourParticleCrossCheckArray[1] = n;
1056  fourParticleCrossCheckArray[2] = -1*m;
1057  fourParticleCrossCheckArray[3] = -1*n;
1058  fourParticleCrossCheckComplex = (CalculateRecursion(4, fourParticleCrossCheckArray))/fourParticleDenominator;
1059  fourParticleCrossCheckValue = fourParticleCrossCheckComplex.Re();
1060 
1061  // Fill the corresponding bin in the TProfile.
1062  iBinFourCrossCheckMiddle = (1.*iBinFourCrossCheck) - 0.5;
1063  fProfileFourParticleCorrelationsCrossCheck->Fill(iBinFourCrossCheckMiddle, fourParticleCrossCheckValue, fourParticleDenominator);
1064  fProfileFourParticleCorrelationsCrossCheck->GetXaxis()->SetBinLabel(iBinFourCrossCheck, Form("(%d,%d)", m, n));
1065 
1066  // Reset of the variables for the next pair of harmonics.
1067  iBinFourCrossCheck++;
1068  fourParticleCrossCheckComplex = TComplex(0., 0.);
1069  fourParticleCrossCheckValue = 0.;
1070  iBinFourCrossCheckMiddle = 0.;
1071  }
1072  }
1073  // Reset the variables for the next event.
1074  iBinFourCrossCheck =1;
1075  }
1076 
1077 // Compute the 6-particle correlations.
1078  Int_t sixZerosArray[6] = {0}; // Zero array for the denominator.
1079  Double_t sixParticleDenominator = CalculateRecursion(6, sixZerosArray).Re(); // 6-particle event weight.
1080 
1081  Int_t sixParticleArray[6] = {0}; // Array to save the harmonics with the format: cos(lphi1 + mphi2 + nphi3 - lphi4 - mphi5 - nphi6).
1082  TComplex sixParticleComplex = TComplex(0., 0.); // Complex value for the 6-particle correlation.
1083  Double_t sixParticleValue = 0.; // Real part of the 6-particle correlation.
1084  Int_t iBinSix = 1; // Index of the corresponding bin in the TProfile.
1085  Int_t iBinSixMiddle = 0.; // Middle of the corresponding bin in the TProfile.
1086  for (Int_t l = 1; l < 5; l++)
1087  {
1088  for (Int_t m = 2; m < 6; m++)
1089  {
1090  if (l >= m) {continue;}
1091  for (Int_t n = 3; n < 7; n++)
1092  {
1093  if ((l >= n) || (m >= n)) {continue;}
1094 
1095  // Compute the correlation.
1096  sixParticleArray[0] = l;
1097  sixParticleArray[1] = m;
1098  sixParticleArray[2] = n;
1099  sixParticleArray[3] = -1*l;
1100  sixParticleArray[4] = -1*m;
1101  sixParticleArray[5] = -1*n;
1102  sixParticleComplex = (CalculateRecursion(6, sixParticleArray))/sixParticleDenominator;
1103  sixParticleValue = sixParticleComplex.Re();
1104 
1105  // Fill the corresponding bin in the TProfile.
1106  iBinSixMiddle = (1.*iBinSix) - 0.5;
1107  fProfileSixParticleCorrelations->Fill(iBinSixMiddle, sixParticleValue, sixParticleDenominator);
1108  fProfileSixParticleCorrelations->GetXaxis()->SetBinLabel(iBinSix, Form("(%d,%d,%d)", l, m, n));
1109 
1110  // Reset of the variables for the next pair of harmonics.
1111  iBinSix++;
1112  sixParticleComplex = TComplex(0., 0.);
1113  sixParticleValue = 0.;
1114  iBinSixMiddle = 0.;
1115  }
1116  }
1117  }
1118 
1119 // Reset the variables for the next event.
1120  twoParticleDenominator = 0.;
1121  fourParticleDenominator = 0.;
1122  sixParticleDenominator = 0.;
1123  iBinFour = 1;
1124  iBinSix = 1;
1125 }
1126 
1127 //======================================================================================//
1129 {
1130 /* Calculate the multi-particle correlators by using the recursion method (an improved faster version) originally developed by Kristjan Gulbrandsen (gulbrand@nbi.dk). */
1131  Int_t nm1 = n-1;
1132  TComplex c(Q(harmonic[nm1], mult));
1133  if (nm1 == 0) return c;
1134  c *= CalculateRecursion(nm1, harmonic);
1135  if (nm1 == skip) return c;
1136 
1137  Int_t multp1 = mult+1;
1138  Int_t nm2 = n-2;
1139  Int_t counter1 = 0;
1140  Int_t hhold = harmonic[counter1];
1141  harmonic[counter1] = harmonic[nm2];
1142  harmonic[nm2] = hhold + harmonic[nm1];
1143  TComplex c2(CalculateRecursion(nm1, harmonic, multp1, nm2));
1144  Int_t counter2 = n-3;
1145 
1146  while (counter2 >= skip) {
1147  harmonic[nm2] = harmonic[counter1];
1148  harmonic[counter1] = hhold;
1149  ++counter1;
1150  hhold = harmonic[counter1];
1151  harmonic[counter1] = harmonic[nm2];
1152  harmonic[nm2] = hhold + harmonic[nm1];
1153  c2 += CalculateRecursion(nm1, harmonic, multp1, counter2);
1154  --counter2;
1155  }
1156  harmonic[nm2] = harmonic[counter1];
1157  harmonic[counter1] = hhold;
1158 
1159  if (mult == 1) return c-c2;
1160  return c-Double_t(mult)*c2;
1161 }
1162 
1163 //======================================================================================//
1164 void AliAnalysisTaskTwoMultiCorrelations::ComputeTwoNestedLoops(long long nParticles, Int_t *harmonic, Double_t aAngles[], Double_t weights[], TProfile *profile, Double_t middleBin)
1165 {
1166 /* Compute the 2-particle correlations using nested loops. */
1167  Double_t twoParticleCosine = 0.; // cos(kphi_1 + lphi_2)).
1168  Double_t twoParticleWeight = 1.; // Total particle weights.
1169 
1170  for (long long m = 0; m < nParticles; m++)
1171  {
1172  for (long long n = 0; n < nParticles; n++)
1173  {
1174  if (m == n) {continue;} // Remove the autocorrelations.
1175 
1176  twoParticleCosine = TMath::Cos(harmonic[0]*aAngles[m] + harmonic[1]*aAngles[n]);
1177  twoParticleWeight = weights[m] + weights[n];
1178  profile->Fill(middleBin, twoParticleCosine, twoParticleWeight);
1179 
1180  // Reset the variables.
1181  twoParticleCosine = 0.;
1182  twoParticleWeight = 1.;
1183  }
1184  }
1185 }
1186 
1187 //======================================================================================//
1188 void AliAnalysisTaskTwoMultiCorrelations::ComputeFourNestedLoops(long long nParticles, Int_t *harmonic, Double_t aAngles[], Double_t weights[], TProfile *profile, Double_t middleBin)
1189 {
1190 /* Compute the 4-particle correlations using nested loops. */
1191  Double_t fourParticleCosine = 0.; // cos(kphi_1 +l phi_2 + mphi_3 + nphi_4).
1192  Double_t fourParticleWeight = 1.; // Total particle weights.
1193 
1194  for (long long k = 0; k < nParticles; k++)
1195  {
1196  for (long long l = 0; l < nParticles; l++)
1197  {
1198  if (k == l) {continue;} // Remove the autocorrelations.
1199  for (long long m = 0; m < nParticles; m++)
1200  {
1201  if ((k == m) || (l == m)) {continue;}
1202  for (long long n = 0; n < nParticles; n++)
1203  {
1204  if ((k == n) || (l == n) || (m == n)) {continue;}
1205 
1206  fourParticleCosine = TMath::Cos(harmonic[0]*aAngles[k] + harmonic[1]*aAngles[l] + harmonic[2]*aAngles[m] + harmonic[3]*aAngles[n]);
1207  fourParticleWeight = weights[k] + weights[l] + weights[m] + weights[n];
1208  profile->Fill(middleBin, fourParticleCosine, fourParticleWeight);
1209 
1210  // Reset the variables.
1211  fourParticleCosine = 0.;
1212  fourParticleWeight = 1.;
1213  }
1214  }
1215  }
1216  }
1217 }
1218 
1219 //######################################################################################//
1220 // Methods called in 'Terminate'.
1221 //======================================================================================//
1222 
1223 
TProfile * fProfileTwoParticleCorrelationsNestedLoops
<6>_{j,k,l,-j,-k,-l} for j = 1..4, k = 2..5 (k > j), l = 3..6 (l > k). (20 bins)
TH1D * fHistoPtAfterSelection
Distribution of the transverse momentum before the track selection.
double Double_t
Definition: External.C:58
TH1I * fHistoInitialNumberOfTracks
Distribution of the centrality of the events.
TH1D * fHistoDCAZAfterSelection
Distribution of the z-coordinate of the DCA before the track selection.
virtual void ComputeFourNestedLoops(long long nParticles, Int_t *harmonic, Double_t aAngles[], Double_t weights[], TProfile *profile, Double_t middleBin)
TH1D * fHistoEtaBeforeSelection
Distribution of the transverse momentum after the full selection.
TH1D * fHistoTPCChiSquareAfterSelection
Distribution of the chi square of the track momentum in the TPC before the track selection.
TH1I * fHistoTPCClustersBeforeSelection
Distribution of the azimuthal angles after the full selection.
TH1D * fHistoPtBeforeSelection
Distribution of the PV z-position after the full selection.
TProfile * fProfileTwoParticleCorrelations
Distribution of the z-coordinate of the DCA after the full selection.
TCanvas * c
Definition: TestFitELoss.C:172
TH1I * fHistoNumberOfTracksBeforeTrackSelection
Distribution of the initial number of tracks.
TH1D * fHistoDCAZBeforeSelection
Distribution of the xy-plane of the DCA after the full selection.
TProfile * fProfileFourParticleCorrelationsNestedLoops
<2>_{j,-j} for j = 1..6 with nested loops. (6 bins)
int Int_t
Definition: External.C:63
Definition: External.C:204
TH1I * fHistoTPCClustersAfterSelection
Distribution of the number of TPC clusters before the track selection.
TH1D * fHistoVertexZBeforeSelection
Distribution of the PV y-position after the full selection.
TH1D * fHistoEtaAfterSelection
Distribution of the pseudorapidity before the track selection.
TH1D * fHistoPhiAfterSelection
Distribution of the azimuthal angles before the track selection.
Definition: External.C:212
virtual void CalculateQvectors(long long numberOfParticles, Double_t angles[], Double_t pWeights[])
TH1D * fHistoDCAXYBeforeSelection
Distribution of the chi square of the track momentum in the TPC after the full selection.
TH1D * fHistoVertexZAfterSelection
Distribution of the initial PV z-position.
virtual void ComputeTwoNestedLoops(long long nParticles, Int_t *harmonic, Double_t aAngles[], Double_t weights[], TProfile *profile, Double_t middleBin)
TH1D * fHistoVertexYAfterSelection
Distribution of the initial PV y-position.
TH1D * fHistoDCAXYAfterSelection
Distribution of the xy-plane of the DCA before the track selection.
TProfile * fProfileFourParticleCorrelationsCrossCheck
<4>_{j,k,-j,-k} for j = 1..6, k = j..6. (21 bins)
TH1D * fHistoPhiBeforeSelection
Distribution of the pseudorapidity after the full selection.
TH1D * fHistoVertexYBeforeSelection
Distribution of the PV x-position after the full selection.
TH1D * fHistoTPCChiSquareBeforeSelection
Distribution of the number of TPC clusters after the full selection.
TH1I * fHistoFinalNumberOfTracks
Distribution of the number of tracks before the track selection.
TProfile * fProfileFourParticleCorrelations
<2>_{j,-j} for j = 1..6. (6 bins)
TProfile * fProfileSixParticleCorrelations
<4>_{j,k,-j,-k} for j = 2..6, k = 1..j-1. (15 bins, cross-check purpose)
TComplex CalculateRecursion(Int_t n, Int_t *harmonic, Int_t mult=1, Int_t skip=0)
virtual void AnalyseAODevent(AliAODEvent *aAODevent)
const char Option_t
Definition: External.C:48
TH1D * fHistoVertexXAfterSelection
Distribution of the initial PV x-position.
bool Bool_t
Definition: External.C:53
virtual void ComputeMultiparticleCorrelations(long long numberOfParticles, Double_t angles[], Double_t pWeights[])
Bool_t ApplyTrackSelection(Double_t momentum, Double_t pseudorapidity, Int_t NclustersInTPC, Double_t TPCchiSquare, Double_t xyDCA, Double_t zDCA)