AliPhysics  9c66e61 (9c66e61)
AliAnalysisTaskTwoMultiCorrelations.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 
16 /*******************************************************************************
17 * Analysis task for anisotropic flow analysis of data taken by ALICE *
18 * with different methods for two- and multiparticle correlations *
19 * *
20 * Author: Cindy Mordasini (cindy.mordasini@cern.ch) *
21 *******************************************************************************/
22 
23 #include "Riostream.h"
24 #include <vector>
26 #include "AliLog.h"
27 #include "AliAODEvent.h"
28 #include "AliAODInputHandler.h"
29 #include "AliAnalysisManager.h"
30 #include "AliMultSelection.h"
31 #include "TFile.h"
32 #include "TComplex.h"
33 #include "TMath.h"
34 
35 using std::cout;
36 using std::endl;
37 
39 
40 // ============================================================================================== //
41 
43  AliAnalysisTaskSE(name),
44  fNparticlesCorrelations(2), // Number of m-particle correlations and harmonics (2-14)
45  fHarmonicOne(2), // Harmonic n_1, default value for 2-p correlation
46  fHarmonicTwo(-2), // Harmonic n_2, default value for 2-p correlation
47  fHarmonicThree(0), // Harmonic n_3
48  fHarmonicFour(0), // Harmonic n_4
49  fHarmonicFive(0), // Harmonic n_5
50  fHarmonicSix(0), // Harmonic n_6
51  fHarmonicSeven(0), // Harmonic n_7
52  fHarmonicEight(0), // Harmonic n_8
53  fHarmonicNine(0), // Harmonic n_9
54  fHarmonicTen(0), // Harmonic n_10
55  fHarmonicEleven(0), // Harmonic n_11
56  fHarmonicTwelve(0), // Harmonic n_12
57  fHarmonicThirteen(0), // Harmonic n_13
58  fHarmonicFourteen(0), // Harmonic n_14
59  fMinCentrality(0.0), // Minimum of centrality
60  fMaxCentrality(100.0), // Maximum of centrality
61  fUseParticleWeights(kFALSE), // Use non-unit particle weights
62  fDoNestedLoops(kFALSE), // Cross-check the results with nested loops
63  fOutputList(NULL), // Main list holding all the output objects
64  fControlOutputList(NULL), // List holding all the control objects
65  fDraftOutputList(NULL), // List holding all the intermediate objects
66  fFinalOutputList(NULL), // List holding all the final results
67  fPtControlHisto(NULL), // Control histogram for the transverse momentum
68  fEtaControlHisto(NULL), // Control histogram for the pseudorapidity
69  fControlPhiHisto(NULL), // Control histogram for the azimuthal angles
70  fCentralityHisto(NULL), // Control histogram for the centrality
71  fMultiplicityDist(NULL), // Control histogram for the multiplicity distribution
72  fCorrelationWithQvectorsProfile(NULL), // m-particle correlation estimated with Q-vectors
73  fCorrelationWithNestedLoopsProfile(NULL), // 2-p correlation estimated with nested loops
74  fCorrelationWithQvectorsSaProfile(NULL) // 2-p correlation estimated with stand-alone Q-vectors
75  //fEstimatedFlowWithQcProfile(NULL) // Anisotropic flow estimated with Q-cumulants
76 {
77 // Constructor of the class
78 
79  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
80 
81  // Creation of a new main list
82  fOutputList = new TList();
83  fOutputList->SetName("outputAnalysis");
84  fOutputList->SetOwner(kTRUE);
85 
86  // Definition of the input and output slots
87  DefineOutput(1, TList::Class());
88 
89  if(useParticleWeights)
90  {
91  // not needed for the time being, maybe insert here the call for the file with particle weights???
92  }
93 
94 } // End of the constructor
95 
96 // ********************************************************************************************** //
97 
100  fNparticlesCorrelations(2), // Number of m-particle correlations and harmonics (2-14)
101  fHarmonicOne(2), // Harmonic n_1
102  fHarmonicTwo(-2), // Harmonic n_2
103  fHarmonicThree(0), // Harmonic n_3
104  fHarmonicFour(0), // Harmonic n_4
105  fHarmonicFive(0), // Harmonic n_5
106  fHarmonicSix(0), // Harmonic n_6
107  fHarmonicSeven(0), // Harmonic n_7
108  fHarmonicEight(0), // Harmonic n_8
109  fHarmonicNine(0), // Harmonic n_9
110  fHarmonicTen(0), // Harmonic n_10
111  fHarmonicEleven(0), // Harmonic n_11
112  fHarmonicTwelve(0), // Harmonic n_12
113  fHarmonicThirteen(0), // Harmonic n_13
114  fHarmonicFourteen(0), // Harmonic n_14
115  fMinCentrality(0.0), // Minimum of centrality
116  fMaxCentrality(100.0), // Maximum of centrality
117  fUseParticleWeights(kFALSE), // Use non-unit particle weights
118  fDoNestedLoops(kFALSE), // Cross-check the results with nested loops
119  fOutputList(NULL), // Main list holding all the output objects
120  fControlOutputList(NULL), // List holding all the control objects
121  fDraftOutputList(NULL), // List holding all the intermediate objects
122  fFinalOutputList(NULL), // List holding all the final results
123  fPtControlHisto(NULL), // Control histogram for the transverse momentum
124  fEtaControlHisto(NULL), // Control histogram for the pseudorapidity
125  fControlPhiHisto(NULL), // Control histogram for the azimuthal angles
126  fCentralityHisto(NULL), // Control histogram for the centrality
127  fMultiplicityDist(NULL), // Control histogram for the multiplicity distribution
128  fCorrelationWithQvectorsProfile(NULL), // m-particle correlation estimated with Q-vectors
129  fCorrelationWithNestedLoopsProfile(NULL), // 2-p correlation estimated with nested loops
130  fCorrelationWithQvectorsSaProfile(NULL) // 2-p correlation estimated with stand-alone Q-vectors
131  //fEstimatedFlowWithQcProfile(NULL) // Anisotropic flow estimated with Q-cumulants
132 {
133 // Dummy constructor of the class
134 
135  AliDebug(2, "AliAnalysisTaskTwoMultiCorrelations::AliAnalysisTaskTwoMultiCorrelations(const char *name, Bool_t useParticleWeights)");
136 
137 } // End of the dummy constructor
138 
139 // ********************************************************************************************** //
140 
142 {
143 // Destructor of the class
144  // Delete the main TList => delete automatically everything holds in it
145 
146  if(fOutputList) delete fOutputList;
147 
148 } // End of the destructor
149 
150 // ============================================================================================== //
151 
153 {
154 // Method called at every worker node to initialise the lists
155 // Organisation of the method
156  // First part of the trick to avoid name clashes
157  // Booking and nesting of all the lists
158  // Booking of all the objects
159  // Second part of the trick to avoid name clashes
160 
161 // First part of the trick to avoid name clashes
162  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
163  TH1::AddDirectory(kFALSE);
164 
165 // Booking and nesting of all the lists
166  this->BookAndNestAllLists();
167 
168 // Booking of all the objects
169  this->BookControlList();
170  this->BookDraftList();
171  this->BookFinalList();
172 
173 // Second part of the trick to avoid name clashes
174  TH1::AddDirectory(oldHistAddStatus);
175 
176  PostData(1,fOutputList);
177 
178 } // End of AliAnalysisTaskTwoMultiCorrelations::UserCreateOutputObjects()
179 
180 // ********************************************************************************************** //
181 
183 {
184 // Method called for each event, contains all the calculations
185  // Note to self: find a way to include non-unit particle weights and select them via fUseParticleWeights
186 // Organisation of the method
187  // Obtention of the pointer to the AOD event
188  // Gestion of the centrality
189  // Loop over the multiplicity of the event
190  // Obtention of the pointer to the current particle
191  // Definition and control histograms of some variables for the cuts
192  // Determination of the number of particles that pass the cuts
193  // Definition of the varibles for the analysis
194  // Do the analysis only if the number of particles after the cuts is higher than the one for the m-particle correlation
195  // Gestion of the centrality
196  // Loop over the number of remaining particles
197  // Obtention of the pointer to the current particle
198  // Filling of the azimuthal angles and particle weights and the control histograms
199  // Computation of the correlations with different possible methods
200  // Release of the allocated memory
201  // PostData
202 
203 //cout << "TEST" << endl;
204 //exit(0);
205 
206 
207 // Obtention of the pointer to the AOD event (from TaskSE)
208  AliAODEvent *currentEvent = dynamic_cast<AliAODEvent*>(InputEvent());
209  if(!currentEvent){return;} // Protection against NULL pointer
210 
211 // Gestion of the centrality
212  AliMultSelection *ams = (AliMultSelection*)currentEvent->FindListObject("MultSelection");
213  if(!ams){return;}
214 
215 // Loop over the multiplicity of the event
216  Int_t nParticles = currentEvent->GetNumberOfTracks(); // Number of particles of the event
217  Int_t nParticlesAfterCuts = 0; // Number of particles that pass the various cuts (pT, eta, ...)
218  Int_t *goodIndices = new Int_t[nParticles](); // Vector where the index of a particle which passes the cut is saved with 1 and if it does not pass, the index corresponds to 0
219 
220  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++)
221  {
222  // Obtention of the pointer to the current particle
223  AliAODTrack *currentParticle = dynamic_cast<AliAODTrack*>(currentEvent->GetTrack(iParticle));
224  if(!currentParticle){continue;} // Protection against NULL pointer
225  if(!currentParticle->TestFilterBit(128)){continue;} // Filter bit 128 denotes TPC-only tracks
226 
227  // Definition and control histograms of some variables for the cuts
228  Double_t pT = currentParticle->Pt(); // Transverse momentum
229  Double_t eta = currentParticle->Eta(); // Pseudorapidity
230 
231  fPtControlHisto->Fill(pT);
232  fEtaControlHisto->Fill(eta);
233 
234  // Determination of the number of particles which pass the cuts and set the flag 1/0 of goodIndices[iParticle]
235  if ( (-0.8 < eta) && (eta < 0.8) && (0.2 < pT) && (pT < 5.0) )
236  {
237  nParticlesAfterCuts++;
238  goodIndices[iParticle] = 1;
239  }
240  else {goodIndices[iParticle] = 0;}
241 
242  } // End of for (Int_t iParticle = 0; iParticle < nParticles; iParticle++)
243 
244 // Filling of the control histogram for the multiplicity distribution
245  fMultiplicityDist->Fill(nParticlesAfterCuts);
246 
247 // Definition of the varibles for the analysis
248  Double_t *phi = new Double_t[nParticlesAfterCuts](); // Azimuthal angles
249  Double_t *particleWeights = new Double_t[nParticlesAfterCuts](); // Particle weights
251  Int_t index = 0; // Index of the "good index", increased when the loop reaches a 1 in goodIndices
252 
253 // Do the analysis only if the number of particles after the cuts is higher than the one for the m-particle correlation
254  // (in order to avoid the division by zero)
255  if (nParticlesAfterCuts >= fNparticlesCorrelations)
256  {
257  // Gestion of the centrality
258  if(ams->GetMultiplicityPercentile("V0M") >= fMinCentrality && ams->GetMultiplicityPercentile("V0M") < fMaxCentrality)
259  {
260  fCentralityHisto->Fill(ams->GetMultiplicityPercentile("V0M"));
261  }
262  else {return;} // this event does not belong to the centrality class specified for this particular analysis
263 
264  // Loop over the number of remaining particles
265  for (Int_t iiParticle = 0; iiParticle < nParticles; iiParticle++)
266  {
267  // Obtention of the pointer to the current particle
268  AliAODTrack *keptParticle = dynamic_cast<AliAODTrack*>(currentEvent->GetTrack(iiParticle));
269  if(!keptParticle){continue;} // Protection against NULL pointer
270  if(!keptParticle->TestFilterBit(128)){continue;} // Filter bit 128 denotes TPC-only tracks
271 
272  // Filling of the azimuthal angles and particle weights and the control histograms only if goodIndices[iiParticle] == 1
273  if (goodIndices[iiParticle] == 1)
274  {
275  phi[index] = keptParticle->Phi();
276  particleWeights[index] = 1.;
277  //if (fUseParticleWeights) {continue;} // Note to self: add the gestion of non-unit weight from external file
278  //else {particleWeights[iiParticle] = 1.;}
279 
280  fControlPhiHisto->Fill(phi[index]);
281  index++;
282  }
283  else {continue;} // End of if ( (-0.8 < eta) && (eta < 0.8) && (0.2 < pT) && (pT < 5.0) )
284  } // End of for (Int_t iiParticle = 0; iiParticle < nParticlesAfterCuts; iiParticle++)
285 
286  // Computation of the correlations with different possible methods
287  // Q-vectors with recursion formula
288  for (Int_t iParticleCorr = 2; iParticleCorr <= fNparticlesCorrelations; iParticleCorr = iParticleCorr + 2)
289  {
290  //ComputeCorrelationsWithQvectors(nParticles, phi, particleWeights, harmonics, fNparticlesCorrelations);
291  ComputeCorrelationsWithQvectors(nParticlesAfterCuts, phi, particleWeights, harmonics, iParticleCorr);
292  } // End of for (Int_t iCorr = 2; iCorr <= fNparticlesCorrelations; iCorr = iCorr + 2)
293 
294  // Nested loops (for cross-checking the results from the Q-vectors)
295  if ((fDoNestedLoops) && (fNparticlesCorrelations == 2)) {ComputeCorrelationsWithTwoNestedLoops(fHarmonicOne, fHarmonicTwo, nParticlesAfterCuts, phi, particleWeights);}
297 
298  // Q-vectors with stand-alone formula pour 2-p correlation (debug for recursion)
299  if (fNparticlesCorrelations == 2) {ComputeCorrelationsWithStandAloneQvectors(fHarmonicOne, 0, nParticlesAfterCuts, phi, particleWeights);}
300 
301  // Release of the allocated memory
302  delete [] goodIndices;
303  delete [] phi;
304  delete [] particleWeights;
305 
306  } // End of if (nParticlesAfterCuts >= fNparticlesCorrelations)
307 
308 // PostData
309  PostData(1,fOutputList);
310 
311 } // End of void AliAnalysisTaskTwoMultiCorrelations::UserExec(Option_t *)
312 
313 //******************************************************************************
314 
316 {
317 // Method called at the end of the execution, once the run over the events is over
318 // Organisation of the method
319  // 1.) Access to the merged output list
320  // 2.) Estimation of anisotropic flow with Q-cumulants
321  // 3.) Creation of the output file and save of the main list in it
322 
323 // 1.) Access to the merged output list
324  fOutputList = (TList*)GetOutputData(1);
325  if(!fOutputList){exit(1);}
326 
327 // 2.) Estimation of anisotropic flow with Q-cumulants (CANNOT BE MERGED)
328  //EstimateFlowWithCumulants(fNparticlesCorrelations);
329 
330 // 3.) Creation of the output file and save of the main list in it
331  TFile *outputFile = new TFile("AnalysisResults.root", "RECREATE");
332  fOutputList->Write(fOutputList->GetName(),TObject::kSingleKey);
333  delete outputFile;
334 
335 } // End of void AliAnalysisTaskTwoMultiCorrelations::Terminate(Option_t *)
336 
337 //==============================================================================
338 
340 {
341 // Method to prepare all the lists with the results in the output file
342 // Organisation of the method
343  // 1.) Booking and nesting lists for the control objects
344  // 2.) Booking and nesting lists for the intermediate objects
345  // 3.) Booking and nesting lists for the final results
346 
347  TString sMethodName = "void AliAnalysisTaskTwoMultiCorrelations::BookAndNestAllLists()";
348  if(!fOutputList){Fatal(sMethodName.Data(),"Main list fOutputList is NULL");}
349 
350 // 1.) Booking and nesting lists for the control objects
351  fControlOutputList = new TList();
352  fControlOutputList->SetName("ControlOutputList");
353  fControlOutputList->SetOwner(kTRUE);
354  fOutputList->Add(fControlOutputList);
355 
356 // 2.) Booking and nesting lists for the intermediate objects
357  fDraftOutputList = new TList();
358  fDraftOutputList->SetName("DraftOutputList");
359  fDraftOutputList->SetOwner(kTRUE);
361 
362 // 3.) Booking and nesting lists for the final results
363  fFinalOutputList = new TList();
364  fFinalOutputList->SetName("FinalOutputList");
365  fFinalOutputList->SetOwner(kTRUE);
367 
368 } // End of AliAnalysisTaskTwoMultiCorrelations::BookAndNestAllLists()
369 
370 //==============================================================================
371 
373 {
374 // Method to prepare the list with the control histograms
375 // Organisation of the method
376  // Control histogram for the transverse momentum
377  // Control histogram for the pseudorapidity
378  // 1.) Control histogram for the azimuthal angles
379 
380 // Control histrogram for the transverse momentum
381  fPtControlHisto = new TH1F("fPtControlHisto", "Transverse momentum distribution", 1000, 0., 20.);
382  fPtControlHisto->SetStats(kTRUE);
383  fPtControlHisto->GetXaxis()->SetTitle("p_{T}");
385 
386 // Control histogram for the pseudorapidity
387  fEtaControlHisto = new TH1F("fEtaControlHisto", "Pseudorapidity distribution", 1000, -1., 1.);
388  fEtaControlHisto->SetStats(kTRUE);
389  fEtaControlHisto->GetXaxis()->SetTitle("eta");
391 
392 // 1.) Control histogram for the azimuthal angles
393  fControlPhiHisto = new TH1F("fControlPhiHisto", "Azimuthal angles distribution", 1000, 0., 6.3);
394  fControlPhiHisto->SetStats(kTRUE);
395  fControlPhiHisto->GetXaxis()->SetTitle("phi");
397 
398 // 2.) Control histogram for the centrality
399  fCentralityHisto = new TH1F("fCentralityHisto", "Centrality distribution", 9, 0., 100.);
400  fCentralityHisto->SetStats(kTRUE);
401  fCentralityHisto->GetXaxis()->SetTitle("Centrality percentile");
403 
404 // 3.) Control TProfile for the average multiplicity
405  fMultiplicityDist = new TH1F("fMultiplicityDist", "Multiplicity distribution", 3000, 0., 3000.);
406  fMultiplicityDist->GetXaxis()->SetTitle("Multiplicity");
408 
409 } // End of void AliAnalysisTaskTwoMultiCorrelations::BookControlList()
410 
411 //******************************************************************************
412 
414 {
415 // Method to prepare the list with the intermediate results from the computation of Q-vectors and nested loops
416 // Organisation of the method
417  // 1.) TProfile from the method of the Q-vectors
418  // 2.) TProfile from the method of the two nested loops
419  // 3.) TProfile from the method of the Q-vectors, stand alone formula
420 
421 // 1.) TProfile from the method of the Q-vectors
422  fCorrelationWithQvectorsProfile = new TProfile("fCorrelationWithQvectorsProfile", "m-particle correlation with Q-vectors", static_cast<int>(fNparticlesCorrelations)/2, 0., static_cast<double>(fNparticlesCorrelations)/2.);
423  fCorrelationWithQvectorsProfile->SetStats(kFALSE);
425  //fCorrelationWithQvectorsProfile->GetXaxis()->SetTitle("m-particle correlation");
427 
428 // 2.) TProfile from the method of nested loops
429  fCorrelationWithNestedLoopsProfile = new TProfile("fCorrelationWithNestedLoopsProfile", "m-particle correlation with nested loops", 1, 0., 1.);
430  fCorrelationWithNestedLoopsProfile->SetStats(kFALSE);
432  //fCorrelationWithNestedLoopsProfile->GetXaxis()->SetTitle("m-particle correlation");
434 
435 // 3.) TProfile from the method of the Q-vectors, stand alone formula
436  fCorrelationWithQvectorsSaProfile = new TProfile("fCorrelationWithQvectorsSaProfile", "m-particle correlation with Q-vectors, stand alone", 1, 0., 1.);
437  fCorrelationWithQvectorsSaProfile->SetStats(kTRUE);
439  //fCorrelationWithQvectorsSaProfile->GetXaxis()->SetTitle("2-particle correlation");
441 
442 } // End of void AliAnalysisTaskTwoMultiCorrelations::BookDraftList()
443 
444 //******************************************************************************
445 
447 {
448 // Method to prepare the list with the final results for the anisotropic flow
449 // Organisation of the method
450  // 1.) TProfile of the anisotropic flow estimated with Q-cumulants
451  // 2.) TProfile of the errors on the anisotropic flow estimated with Q-cumulants
452 
453 // 1.) TProfile of the anisotropic flow estimated with Q-cumulants
454  /*fEstimatedFlowWithQcProfile = new TProfile("fEstimatedFlowWithQcProfile", "Flow v_n estimated with Q-cumulants", 1, 0., 1.);
455  fEstimatedFlowWithQcProfile->SetStats(kTRUE);
456  //fEstimatedFlowWithQcProfile->Sumw2();
457  fEstimatedFlowWithQcProfile->GetXaxis()->SetTitle("Estimated v_n");
458  fFinalOutputList->Add(fEstimatedFlowWithQcProfile);*/
459 
460 // 2.) TProfile of the errors on the anisotropic flow estimated with Q-cumulants
461  /*fEstimatedFlowWithQcErrorsProfile = new TProfile("fEstimatedFlowWithQcErrorsProfile", "Error on v_n estimated with QC", 1, 0., 1.);
462  fEstimatedFlowWithQcErrorsProfile->SetStats(kTRUE);
463  fEstimatedFlowWithQcErrorsProfile->GetXaxis()->SetTitle("Error on estimated v_n");
464  fFinalOutputList->Add(fEstimatedFlowWithQcErrorsProfile);*/
465 
466 } // End of void AliAnalysisTaskTwoMultiCorrelations::BookFinalList()
467 
468 //==============================================================================
469 
470 TComplex AliAnalysisTaskTwoMultiCorrelations::CalculateQvector(Int_t n, Int_t p, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
471 {
472 // Method calculating the "general" definition of the Q-vector Q_(n,p) for arbitrary (n,p)
473 // Organisation of the method
474  // 1.) Declaration of the local variables: Q_(n,p)
475  // 2.) Computation of the Q-vector
476  // 3.) Application of the property Q_(-n,p) = Q_(n,p)* and return of the result
477 
478 // 1.) Declaration of the local variables: Q_(n,p)
479  TComplex qVectorNp = TComplex(0,0); // Q-vector Q_(n,p) for the current event
480  Double_t pWeightPowerP = 0.; // Particle weight to the power p
481 
482 // 2.) Computation of the Q-vector
483  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++)
484  {
485  pWeightPowerP = pow(particleWeight[iParticle], p);
486  qVectorNp += pWeightPowerP * TComplex::Exp((static_cast<double>(n))*(TComplex::I())*phi[iParticle]);
487  } // End of for (Int_t iParticle = 0; iParticle < nParticles; iParticle++)
488 
489 // 3.) Application of the property Q_(-n,p) = Q_(n,p)* and return of the result
490  //if (n < 0) {return TComplex::Conjugate(qVectorNp);}
491  if (n < 0) {CalculateQvector(-n, p, nParticles, phi, particleWeight);}
492  return qVectorNp;
493 
494 } // End of TComplex AliAnalysisTaskTwoMultiCorrelations::CalculateQvector(Int_t n, Int_t p, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
495 
496 //******************************************************************************
497 
498 TComplex AliAnalysisTaskTwoMultiCorrelations::CalculateRecursionWithQvectors(Int_t nParticles, Double_t phi[], Double_t particleWeight[], Int_t nCorr, Int_t harmonics[], Int_t p, Int_t skip)
499 {
500 // Method calculating the recursion with the Q-vectors for the numerator of the m-particle correlation according to the generic framework
501 // Inspired by the method originally developped by Kristjan Gulbrandsen (gulbrand@nbi.dk)
503 // Organisation of the method
504  // 1.) Declaration of the local variables
505  // 2.) Stop conditions of the recursion
506  // 3.) Computation of the recursion
507  // 4.) Return of the result after the recursion
508 
509 // 1.) Declaration of the local variables
510  Int_t nMinusOne = 0; // Harmonic (nCorr-1)
511  TComplex stopQvector = TComplex(0,0);// Q-vector used to stop the recursion
512  Int_t pPlusOne = 0; // (p + 1)
513  Int_t nMinusTwo = 0; // Harmonic (n-2)
514  Int_t counterOne = 0; // First counter for the intermediate indices
515  Int_t hHold = 0; // Temporary harmonic
516  TComplex tempQvector = TComplex(0,0);// Temporary Q-vector
517  Int_t counterTwo = 0; // Second counter for the intermediate indices
518 
519 // 2.) Stop conditions of the recursion
520  nMinusOne = nCorr-1;
521  stopQvector = CalculateQvector(harmonics[nMinusOne], p, nParticles, phi, particleWeight);
522 
523  if (nMinusOne == 0) {return stopQvector;}
524  stopQvector *= CalculateRecursionWithQvectors(nParticles, phi, particleWeight, nMinusOne, harmonics);
525  if (nMinusOne == skip) {return stopQvector;}
526 
527 // 3.) Computation of the recursion
528  pPlusOne = p+1;
529  nMinusTwo = nCorr-2;
530  hHold = harmonics[counterOne];
531  harmonics[counterOne] = harmonics[nMinusTwo];
532  harmonics[nMinusTwo] = hHold + harmonics[nMinusOne];
533  tempQvector = CalculateRecursionWithQvectors(nParticles, phi, particleWeight, nMinusOne, harmonics, pPlusOne, nMinusTwo);
534  counterTwo = nCorr-3;
535 
536  while (counterTwo >= skip)
537  {
538  harmonics[nMinusTwo] = harmonics[counterOne];
539  harmonics[counterOne] = hHold;
540  ++counterOne;
541 
542  hHold = harmonics[counterOne];
543  harmonics[counterOne] = harmonics[nMinusTwo];
544  harmonics[nMinusTwo] = hHold + harmonics[nMinusOne];
545  tempQvector += CalculateRecursionWithQvectors(nParticles, phi, particleWeight, nMinusOne, harmonics, pPlusOne, counterTwo);
546  --counterTwo;
547 
548  } // End of the while (counterTwo >= skip)
549 
550  harmonics[nMinusTwo] = harmonics[counterOne];
551  harmonics[counterOne] = hHold;
552 
553 // 4.) Return of the result after the recursion
554  if (p == 1) {return stopQvector - tempQvector;}
555  return stopQvector - (Double_t(p)*tempQvector);
556 
557 } // End of void AliAnalysisTaskTwoMultiCorrelations::CalculateRecursionWithQvectors(Int_t nParticles, Double_t *phi[], Double_t *particleWeight[], Int_t n, Double_t *harmonics[], Int_t p = 1, Int_t skip = 0)
558 
559 //******************************************************************************
560 
561 void AliAnalysisTaskTwoMultiCorrelations::ComputeCorrelationsWithQvectors(Int_t nParticles, Double_t phi[], Double_t particleWeight[], Int_t harmonics[], Int_t nCorr)
562 {
563 // Method to compute the m-particle correlation using the numerator for the generic framework computed with the recursion method
564 // Organisation of the method
565  // 1.) Declaration of the local variables
566  // 2.) Filling of the array for the numerator
567  // 3.) Computation of the m-particle correlation
568 
569 // 1.) Declaration of the local variables
570  std::vector<Int_t> numeratorHarmonics(nCorr); // Harmonics with values from harmonics array
571  std::vector<Int_t> denominatorHarmonics(nCorr); // Harmonics with value 0 for the denominator
572  memset(numeratorHarmonics.data(), 0, sizeof(Int_t) * numeratorHarmonics.size());
573  memset(denominatorHarmonics.data(), 0, sizeof(Int_t) * denominatorHarmonics.size());
574 
575  Double_t denominator = 0.; // Denominator of <m>_(n1...nm) and event weight for the TProfile
576  TComplex mParticleCorrelation = TComplex(0,0); // m-particle correlation
577  TComplex eventWeight = TComplex(0,0);// Event weight
578 
579  Double_t middleBin = 0.; // Compute the middle value of the bin corresponding to <<m>> in the TProfile
580 
581 // 2.) Filling of the array for the numerator
582  for (Int_t iCorr = 0; iCorr < nCorr; iCorr++)
583  {
584  numeratorHarmonics[iCorr] = harmonics[iCorr];
585  denominatorHarmonics[iCorr] = 0;
586  } // End of for (Int_t iCorr = 0; iCorr < nCorr; iCorr++)
587 
588 // 3.) Computation of the m-particle correlation
589  denominator = (CalculateRecursionWithQvectors(nParticles, phi, particleWeight, nCorr, denominatorHarmonics.data())).Re();
590  cout << Form("Denominator: %f", denominator) << endl;
591  mParticleCorrelation = CalculateRecursionWithQvectors(nParticles, phi, particleWeight, nCorr, numeratorHarmonics.data())/denominator;
592  //eventWeight = (CalculateRecursionWithQvectors(nParticles, phi, particleWeight, nCorr, denominatorHarmonics)).Re();
593 
594 // 4.) Filling of the TProfile
595  middleBin = (static_cast<double>(nCorr) - 1.)/2.;
596  fCorrelationWithQvectorsProfile->Fill(middleBin, mParticleCorrelation.Re(), denominator);
597 
598 } // End of void AliAnalysisTaskTwoMultiCorrelations::ComputeCorrelationsWithQvectors(Int_t n, Int_t nParticles, Double_t *phi[], Double_t *particleWeight[], Int_t *harmonics[], Int_t nCorr)
599 
600 //******************************************************************************
601 
603 {
604 // Method to compute the 2-particle correlation with two nested loops (for cross-checking results)
605 // Organisation of the method
606  // 1.) Declaration of the local variables
607  // 2.) Computation of the 2-p correlation, <2>
608  // 3.) Reset of the local variables before changing the event
609 
610 // 1.) Declaration of the local variables
611  Double_t twoParticleCorrelation = 0.; // Single-event 2-p correlation, <2>_(n1...nm)
612  Double_t totalParticleWeight = 0.; // Particle weight for the 2-p correlation
613 
614 // 2.) Computation of the 2-p correlation, <2>
615  for (Int_t k = 0; k < nParticles; k++)
616  {
617  for (Int_t l = 0; l < nParticles; l++)
618  {
619  // Removal of the autocorrelations k == l
620  if (k == l) continue;
621 
622  // Computation of <2> for a pair of particles
623  twoParticleCorrelation = cos(nOne*phi[k] + nTwo*phi[l]);
624  totalParticleWeight = particleWeight[k] * particleWeight[l];
625 
626  // Filling of the TProfile
627  fCorrelationWithNestedLoopsProfile->Fill(0.5, twoParticleCorrelation, totalParticleWeight);
628 
629  } // End of the loop over the second particle of the pair
630  } // End of the loop over the first particle of the pair
631 
632 // 3.) Reset of the local variables before changing the event
633  twoParticleCorrelation = 0.;
634  totalParticleWeight = 0.;
635 
636 } // End of void AliAnalysisTaskTwoMultiCorrelations::ComputeCorrelationsWithTwoNestedLoops(Int_t nOne, Int nTwo, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
637 
638 //******************************************************************************
639 
641 {
642 // Method to compute the 4-particle correlation with four nested loops (for cross-checking results)
643 // Organisation of the method
644  // 1.) Declaration of the local variables
645  // 2.) Computation of the 4-p correlation, <4>
646  // 3.) Reset of the local variables before changing the event
647 
648 // 1.) Declaration of the local variables
649  Double_t fourParticleCorrelation = 0.; // Single-event 4-p correlation, <4>_(n1,... nm)
650  Double_t totalParticleWeight = 0.;
651 
652 // 2.) Computation of the 4-p correlation, <4>
653  for (Int_t k = 0; k < nParticles; k++)
654  {
655  for (Int_t l = 0; l < nParticles; l++)
656  {
657  // Removal of the autocorrelations k == l
658  if (k == l) continue;
659 
660  for (Int_t m = 0; m < nParticles; m++)
661  {
662  // Removal of the autocorrelations k == m, l == m
663  if ((k == m) || (l == m)) continue;
664 
665  for (Int_t n = 0; n < nParticles; n++)
666  {
667  // Removal of the autocorrelations k == n, l == n, m == n
668  if ((k == n) || (l == n) || (m == n)) continue;
669 
670  // Computation of <4> for a quadruplet of particles
671  fourParticleCorrelation = cos(nOne*phi[k] + nTwo*phi[l] + nThree*phi[m] + nFour*phi[n]);
672  totalParticleWeight = particleWeight[k]*particleWeight[l]*particleWeight[m]*particleWeight[n];
673 
674  // Filling of the TProfile
675  fCorrelationWithNestedLoopsProfile->Fill(0.5, fourParticleCorrelation, totalParticleWeight);
676 
677  } // End of the loop over the fourth particle of the quadruplet
678  } // End of the loop over the third particle of the quadruplet
679  } // End of the loop over the second particle of the quadruplet
680  } // End of the loop over the first particle of the quadruplet
681 
682 // 3.) Reset of the local variables before changing the event
683  fourParticleCorrelation = 0.;
684  totalParticleWeight = 0.;
685 
686 } // End of void AliAnalysisTaskTwoMultiCorrelations::ComputeCorrelationsWithFourNestedLoops(Int_t nOne, Int_t nTwo, Int_t nThree, Int_t nFour, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
687 
688 //******************************************************************************
689 
691 {
692 // Method to compute the 2-p correlation with Q-vectors, stand alone formula (second cross-check for the recursion method
693 // Organisation of the method
694  // 1.) Declaration of the local variables
695  // 2.) Computation of the 2-p correlation
696  // 3.) Filling of the TProfile
697  // 4.) Security: reset to 0 of the variables
698 
699 // 1.) Declaration of the local variables
700  Double_t twoParticleCorrelation = 0.;
701  Int_t twoParticleWeight = 0;
702 
703 // 2.) Computation of the 2-p correlation
704  twoParticleCorrelation = (1./(2.*TMath::Binomial(nParticles, 2))) * (CalculateQvector(n, p, nParticles, phi, particleWeight).Rho2() - nParticles);
705  twoParticleWeight = nParticles * (nParticles - 1);
706 
707 // 3.) Filling of the TProfile
708  fCorrelationWithQvectorsSaProfile->Fill(0.5, twoParticleCorrelation, twoParticleWeight);
709 
710 // 4.) Security: reset of the variables to 0.
711  twoParticleCorrelation = 0.;
712  twoParticleWeight = 0;
713 
714 } // End of void AliAnalysisTaskTwoMultiCorrelations::ComputeCorrelationsWithStandAloneQvectors()
715 
716 //******************************************************************************
717 /*
718 void AliAnalysisTaskTwoMultiCorrelations::EstimateFlowWithCumulants(Int_t maxMcorr)
719 {
720 // Offline method to compute the estimation of v_n{np} using the m-p correlations
721 // obtained with the Q-vectors
722 // Organisation of the method
723  // 1.) Declaration of the local variables
724  // 2.) Filling of the values of the correlations from the TProfile
725  // *.) Computation of the Q-cumulant and the anisotropic flow for all possible cases of maxMcorr
726  // *.) Security: reset of the variables to 0.
727 
728 // 1.) Declaration of the local variables
729  Double_t cNforMaxMcorr = 0.; // Q-cumulant c_n{MaxMcorr}
730  Double_t vNforMaxMcorr = 0.; // Estimated anisotropic flow v_n{MaxMcorr}
731 
732  std::vector<Double_t> mParticleCorrelation(static_cast<int>(maxMcorr/2)); // Vector with the m-particle correlations: <<2>>, <<4>>, <<6>>, <<8>>
733  memset(mParticleCorrelation.data(), 0, sizeof(Double_t) * mParticleCorrelation.size());
734 
735 // 2.) Filling of the values of the correlations from the TProfile
736  for (Int_t i = 0; i < static_cast<int>(maxMcorr/2); i++)
737  {
738  mParticleCorrelation[i] = fCorrelationWithQvectorsProfile->GetBinContent(i+1);
739  } // End of for (Int_t i = 0; i < static_cast<int>(maxMcorr/2); i++)
740 
741 // *.) Computation of the Q-cumulant and the anisotropic flow for all possible cases of maxMcorr
742  switch(maxMcorr)
743  {
744  case 2 : // Computation of c_n{2} and v_n{2}
745  cNforMaxMcorr = mParticleCorrelation[0];
746  vNforMaxMcorr = pow(cNforMaxMcorr, 0.5);
747  break;
748  case 4 :
749  cNforMaxMcorr = mParticleCorrelation[1] - 2.*pow(mParticleCorrelation[0], 2.);
750  vNforMaxMcorr = pow((-1.*cNforMaxMcorr), 0.25);
751  break;
752  case 6 :
753  cNforMaxMcorr = mParticleCorrelation[2] - 9.*mParticleCorrelation[0]*mParticleCorrelation[1] + 12.*pow(mParticleCorrelation[0], 3.);
754  vNforMaxMcorr = pow((cNforMaxMcorr/4.), (1./6.));
755  break;
756  case 8 :
757  cNforMaxMcorr = mParticleCorrelation[3] - 16.*mParticleCorrelation[0]*mParticleCorrelation[2] - 18.*pow(mParticleCorrelation[1], 2.) + 144.*mParticleCorrelation[1]*pow(mParticleCorrelation[2], 2.) - 144.*pow(mParticleCorrelation[0], 4.);
758  vNforMaxMcorr = pow((-1.*cNforMaxMcorr/33.), 0.125);
759  break;
760  } // End of switch
761 
762 // *.) Filling of the TProfile
763  fEstimatedFlowWithQcProfile->Fill(0.5, vNforMaxMcorr);
764  cout << Form("cN: %f", cNforMaxMcorr) << endl;
765 
766 // *.) Security_ reset of the variables to 0.
767  cNforMaxMcorr = 0.;
768  vNforMaxMcorr = 0.;
769  for (Int_t i = 0; i < static_cast<int>(maxMcorr/2); i++)
770  {
771  mParticleCorrelation[i] = 0.;
772  } // End of for (Int_t i = 0; i < static_cast<int>(maxMcorr/2); i++)
773 
774  //return vNforMaxMcorr;
775 
776 }*/ // End of Double_t AliAnalysisTaskTwoMultiCorrelations::ComputeCorrelationsWithStandAloneQvectors()
double Double_t
Definition: External.C:58
virtual void ComputeCorrelationsWithStandAloneQvectors(Int_t n, Int_t p, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
TH1F * fControlPhiHisto
Control histogram for the pseudorapidity.
virtual void ComputeCorrelationsWithQvectors(Int_t nParticles, Double_t phi[], Double_t particleWeight[], Int_t harmonics[], Int_t nCorr)
TProfile * fCorrelationWithQvectorsSaProfile
m-p correlation estimated with nested loops
virtual void ComputeCorrelationsWithTwoNestedLoops(Int_t nOne, Int_t nTwo, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
int Int_t
Definition: External.C:63
TProfile * fCorrelationWithQvectorsProfile
Control histogram for the multiplicity distribution.
virtual void ComputeCorrelationsWithFourNestedLoops(Int_t nOne, Int_t nTwo, Int_t nThree, Int_t nFour, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
TH1F * fMultiplicityDist
Control histogram for the centrality.
TH1F * fEtaControlHisto
Control histogram for the transverse momentum.
TProfile * fCorrelationWithNestedLoopsProfile
m-p correlation estimated with Q-vectors
TComplex CalculateQvector(Int_t n, Int_t p, Int_t nParticles, Double_t phi[], Double_t particleWeight[])
TH1F * fCentralityHisto
Control histogram for the azimuthal angles.
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
TComplex CalculateRecursionWithQvectors(Int_t nParticles, Double_t phi[], Double_t particleWeight[], Int_t nCorr, Int_t harmonics[], Int_t p=1, Int_t skip=0)