AliPhysics  c7b8e89 (c7b8e89)
AliAnalysisTaskStudentsML.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 * template class for student projects *
18 * Marcel Lesch (marcel.lesch@cern.ch) *
19 **********************************************/
20 
21 #include "Riostream.h"
23 #include "AliLog.h"
24 #include "AliAODEvent.h"
25 #include "AliAODInputHandler.h"
26 #include "AliAnalysisManager.h"
27 #include "AliMultSelection.h"
28 #include "TRandom.h"
29 #include "TComplex.h"
30 #include "TProfile.h"
31 #include "TStyle.h"
32 #include "TFile.h"
33 #include "TH1F.h"
34 #include <TArrayF.h>
35 #include <vector>
36 #include "TMath.h"
37 #include "TF1.h"
38 
39 
40 
41 using std::cout;
42 using std::endl;
43 
45 
46 //================================================================================================================
47 
48 AliAnalysisTaskStudentsML::AliAnalysisTaskStudentsML(const char *name, Bool_t useParticleWeights):
49  AliAnalysisTaskSE(name),
50  fHistList(NULL),
51  // Control histograms:
52  fControlHistogramsList(NULL),
53  fPtHist(NULL),
54  fNbins(1000),
55  fMinBin(0.),
56  fMaxBin(10.),
57  fPhiHist(NULL),
58  fEtaHist(NULL),
59  fMultiHisto(NULL),
60  fMaxCorrelator(8),
61  bUseWeights(kFALSE),
62  kNumber(6), //number of correlation
63  kh1(2), kh2(2), kh3(-1), kh4(-1), kh5(-1), kh6(-1), kh7(0), kh8(0), //harmonics
64  ka1(2), ka2(-2), ka3(1), ka4(-1), ka5(1), ka6(-1), ka7(0), ka8(0), //second set of harmonics
65  kSum((kh1<0?-1*kh1:kh1)+(kh2<0?-1*kh2:kh2)+(kh3<0?-1*kh3:kh3)+(kh4<0?-1*kh4:kh4)
66  + (kh5<0?-1*kh5:kh5)+(kh6<0?-1*kh6:kh6)+(kh7<0?-1*kh7:kh7)+(kh8<0?-1*kh8:kh8)), // We will not go beyond 8-p correlations
67  kMaxHarmonic(kSum+1),
68  kMaxPower(fMaxCorrelator+1),
69  fParticles(0),
70  fCentral(0),
71  fMinCentrality(0.),
72  fMaxCentrality(100.),
73  fAngles(NULL),
74  fWeights(NULL),
75  fBin(NULL),
76  func1(NULL),
77  fCentrality(NULL),
78  fCentralitySecond(NULL),
79  fCounterHistogram(NULL),
80  // Final results:
81  fFinalResultsList(NULL)
82  {
83  // Constructor.
84 
85  AliDebug(2,"AliAnalysisTaskStudentsML::AliAnalysisTaskStudentsML(const char *name, Bool_t useParticleWeights)");
86 
87  // Base list:
88  fHistList = new TList();
89  fHistList->SetName("outputStudentAnalysis");
90  fHistList->SetOwner(kTRUE);
91 
92  // Initialize all arrays:
93  this->InitializeArrays();
94 
95  // Define input and output slots here
96  // Input slot #0 works with an AliFlowEventSimple
97  //DefineInput(0, AliFlowEventSimple::Class());
98  // Input slot #1 is needed for the weights input file:
99  //if(useParticleWeights)
100  //{
101  // DefineInput(1, TList::Class());
102  //}
103  // Output slot #0 is reserved
104  // Output slot #1 writes into a TList container
105 
106  DefineOutput(1, TList::Class());
107 
108  if(useParticleWeights)
109  {
110  // not needed for the time being
111  }
112 
113 } // AliAnalysisTaskStudentsML::AliAnalysisTaskStudentsML(const char *name, Bool_t useParticleWeights):
114 
115 //================================================================================================================
116 
119  fHistList(NULL),
120  // Control histograms:
121  fControlHistogramsList(NULL),
122  fPtHist(NULL),
123  fNbins(1000),
124  fMinBin(0.),
125  fMaxBin(10.),
126  fPhiHist(NULL),
127  fEtaHist(NULL),
128  fMultiHisto(NULL),
129  fMaxCorrelator(8),
130  bUseWeights(kFALSE),
131  kNumber(6), //number of correlation
132  kh1(2), kh2(2), kh3(-1), kh4(-1), kh5(-1), kh6(-1), kh7(0), kh8(0), //harmonics
133  ka1(2), ka2(-2), ka3(1), ka4(-1), ka5(1), ka6(-1), ka7(0), ka8(0), //second set of harmonics
134  kSum((kh1<0?-1*kh1:kh1)+(kh2<0?-1*kh2:kh2)+(kh3<0?-1*kh3:kh3)+(kh4<0?-1*kh4:kh4)
135  + (kh5<0?-1*kh5:kh5)+(kh6<0?-1*kh6:kh6)+(kh7<0?-1*kh7:kh7)+(kh8<0?-1*kh8:kh8)), // We will not go beyond 8-p correlations
136  kMaxHarmonic(kSum+1),
137  kMaxPower(fMaxCorrelator+1),
138  fParticles(0),
139  fCentral(0.),
140  fMinCentrality(0.),
141  fMaxCentrality(100.),
142  fAngles(NULL),
143  fWeights(NULL),
144  fBin(NULL),
145  func1(NULL),
146  // Final results:
147  fCentrality(NULL),
148  fCentralitySecond(NULL),
149  fCounterHistogram(NULL),
150  fFinalResultsList(NULL)
151 {
152  // Dummy constructor.
153 
154  AliDebug(2,"AliAnalysisTaskStudentsML::AliAnalysisTaskStudentsML()");
155 
156 } // AliAnalysisTaskStudentsML::AliAnalysisTaskStudentsML():
157 
158 //================================================================================================================
159 
161 {
162  // Destructor.
163 
164  if(fHistList) delete fHistList;
165 
166 } // AliAnalysisTaskStudentsML::~AliAnalysisTaskStudentsML()
167 
168 //================================================================================================================
169 
171 {
172  // Called at every worker node to initialize.
173 
174  // a) Trick to avoid name clashes, part 1;
175  // b) Book and nest all lists;
176  // c) Book all objects;
177  // *) Trick to avoid name clashes, part 2.
178  // d)
179 
180  // a) Trick to avoid name clashes, part 1:
181  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
182  TH1::AddDirectory(kFALSE);
183 
184  // b) Book and nest all lists:
185  this->BookAndNestAllLists();
186 
187  // c) Book all objects:
188  this->BookControlHistograms();
190 
191  // *) Trick to avoid name clashes, part 2:
192  TH1::AddDirectory(oldHistAddStatus);
193 
194  PostData(1,fHistList);
195 
196 } // void AliAnalysisTaskStudentsML::UserCreateOutputObjects()
197 
198 //================================================================================================================
199 
201 {
202  // Main loop (called for each event).
203 
204  // a.0) Get pointer to AOD event;
205  // a.1) Centrality;
206  // b) Start analysis over AODs;
207  // c) Reset event-by-event objects;
208  // d) PostData.
209 
210  fCounterHistogram->Fill(0.5); // counter hist 1st bin
211 
212  // a) Get pointer to AOD event:
213  AliAODEvent *aAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
214  if(!aAOD){return;}
215 
216  fCounterHistogram->Fill(1.5); // counter hist 2nd bin
217 
218  //a.1) Centrality;
219 
220  AliMultSelection *ams = (AliMultSelection*)aAOD->FindListObject("MultSelection");
221  if(!ams){return;}
222 
223  fCounterHistogram->Fill(2.5); // counter hist 3rd bin
224 
225  //func1->SetParameter(2,gRandom->Uniform(0.,TMath::TwoPi())); //*** for testing. for each event psi_2 is different
226 
227 
228  // b) Start analysis over AODs:
229  Int_t nTracks = aAOD->GetNumberOfTracks(); // number of all tracks in current event
230 
231  //fParticles=gRandom->Uniform(50.0,500.0); /*** for testing.
232 
233  Int_t nCounter=0;
234 
235 
236  //first loop: =======================
237 
238  for(Int_t iTrack=0;iTrack<nTracks;iTrack++) // starting a loop over all tracks
239  {
240  AliAODTrack *aTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack)); // getting a pointer to a track
241  if(!aTrack){continue;} // protection against NULL pointers
242  if(!aTrack->TestFilterBit(128)){continue;} // filter bit 128 denotes TPC-only tracks, use only them for the analysis
243 
244  // example variables for each track:
245  Double_t px = aTrack->Px(); // x-component of momenta
246  Double_t py = aTrack->Py(); // y-component of momenta
247  Double_t pz = aTrack->Pz(); // z-component of momenta
248  Double_t e = aTrack->E(); // energy
249  Double_t phi = aTrack->Phi(); // azimuthal angle
250  Double_t eta = aTrack->Eta(); // pseudorapidity
251  Double_t charge = aTrack->Charge(); // charge
252  Double_t pt = aTrack->Pt(); // Pt (transverse momentum)
253 
254 
255  // apply some cuts: e.g. take for the analysis only particles in -0.8 < eta < 0.8, and 0.2 < pT < 5.0
256  // ... implementation of particle cuts ...
257  if( (-0.8 < eta) && (eta < 0.8) && (0.2 < pt) && (pt < 5.0) )
258  {
259  // Fill some control histograms with the particles which passed cuts:
260  fPtHist->Fill(pt);
261  fPhiHist->Fill(phi);
262  fEtaHist->Fill(eta);
263 
264  fParticles += 1; //one more particle passed the cuts
265 
266 
267  // Do some analysis only with the particles which passed the cuts
268  // ... your analysis code ...
269 
270 
271 
272  } // if( (-0.8 < eta) && (eta < 0.8) && (0.2 < pt) && (pt < 5.0) )
273  } // for(Int_t iTrack=0;iTrack<nTracks;iTrack++) // starting a loop over all tracks
274 
275  fAngles = new TArrayD(fParticles);
276 
277  //second loop: =======================
278 
279  for(Int_t iTrack=0;iTrack<nTracks;iTrack++) // starting a loop over all tracks
280  {
281  AliAODTrack *aTrack = dynamic_cast<AliAODTrack*>(aAOD->GetTrack(iTrack)); // getting a pointer to a track
282  if(!aTrack){continue;} // protection against NULL pointers
283  if(!aTrack->TestFilterBit(128)){continue;} // filter bit 128 denotes TPC-only tracks, use only them for the analysis
284 
285  // example variables for each track:
286  Double_t px = aTrack->Px(); // x-component of momenta
287  Double_t py = aTrack->Py(); // y-component of momenta
288  Double_t pz = aTrack->Pz(); // z-component of momenta
289  Double_t e = aTrack->E(); // energy
290  Double_t phi = aTrack->Phi(); // azimuthal angle
291  Double_t eta = aTrack->Eta(); // pseudorapidity
292  Double_t charge = aTrack->Charge(); // charge
293  Double_t pt = aTrack->Pt(); // Pt (transverse momentum)
294 
295 
296  // apply some cuts: e.g. take for the analysis only particles in -0.8 < eta < 0.8, and 0.2 < pT < 5.0
297  // ... implementation of particle cuts ...
298  if( (-0.8 < eta) && (eta < 0.8) && (0.2 < pt) && (pt < 5.0) )
299  {
300 
301  fAngles->AddAt(phi,nCounter);
302  nCounter += 1;
303 
304  } // if( (-0.8 < eta) && (eta < 0.8) && (0.2 < pt) && (pt < 5.0) )
305  } // for(Int_t iTrack=0;iTrack<nTracks;iTrack++) // starting a loop over all tracks
306 
307 
308  //b.1) analysis
309  if(fParticles>0){fMultiHisto->Fill(fParticles);}
310 
311 
312  if(fParticles>=8) //do the correlation only if there are more than 8 particles in the event
313  {
314 
315  if(ams->GetMultiplicityPercentile("V0M") >= fMinCentrality && ams->GetMultiplicityPercentile("V0M") < fMaxCentrality)
316  {
317  if(0. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 5. ){fCentral=0.5;}
318  if(5. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 10. ){fCentral=1.5;}
319  if(10. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 20. ){fCentral=2.5;}
320  if(20. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 30. ){fCentral=3.5;}
321  if(30. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 40. ){fCentral=4.5;}
322  if(40. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 50. ){fCentral=5.5;}
323  if(50. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 60. ){fCentral=6.5;}
324  if(60. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 70. ){fCentral=7.5;}
325  if(70. <= (ams->GetMultiplicityPercentile("V0M")) && (ams->GetMultiplicityPercentile("V0M")) < 80. ){fCentral=8.5;}
326  }
327  else
328  {
329  return; // this event do not belong to the centrality class specified for this particular analysis
330  }
331 
332 
333  Correlation(); //do the correlation
334  fCentrality->Fill(fCentral,fRecursion[0][kNumber-2]->GetBinContent(1),fRecursion[0][kNumber-2]->GetBinContent(2)); //safe output
335  fCentralitySecond->Fill(fCentral,fRecursionSecond[0][kNumber-2]->GetBinContent(1),fRecursionSecond[0][kNumber-2]->GetBinContent(2));
336 
337  fRecursion[0][kNumber-2]->Reset(); //Reset
338  fRecursion[1][kNumber-2]->Reset(); //Reset
339 
340  fRecursionSecond[0][kNumber-2]->Reset(); //Reset
341  fRecursionSecond[1][kNumber-2]->Reset(); //Reset
342 
343 
344 
345 
346  } //if(fParticles>=8)
347 
348 
349 
350 
351 
352  // c) Reset event-by-event objects:
353  fParticles=0;
354  fCentral=0.;
355  nCounter=0;
356  delete fAngles;
357  fAngles=NULL;
358 
359 
360 
361  // d) PostData:
362  PostData(1,fHistList);
363 
364 } // void AliAnalysisTaskStudentsML::UserExec(Option_t *)
365 
366 //================================================================================================================
367 
369 {
370  // Accessing the merged output list.
371 
372  fHistList = (TList*)GetOutputData(1);
373  if(!fHistList){exit(1);}
374 
375  // Do some calculation in offline mode here:
376  // ...
377 
378  TFile *f = new TFile("AnalysisResults.root","RECREATE");
379  fHistList->Write(fHistList->GetName(),TObject::kSingleKey);
380 
381  delete f;
382 
383 } // end of void AliAnalysisTaskStudentsML::Terminate(Option_t *)
384 
385 //================================================================================================================
386 
388 {
389  // Initialize all data members which are arrays in this method.
390 
391  for(Int_t cs=0;cs<2;cs++)
392  {
393  for(Int_t c=0;c<fMaxCorrelator;c++)
394  {
395 
396  fRecursion[cs][c] = NULL;
397  fRecursionSecond[cs][c] = NULL;
398 
399  }
400  } //for(Int_t cs=0;cs<2;cs++)
401 
402  for(Int_t js=0;js<17;js++)
403  {
404  for(Int_t j=0;j<9;j++)
405  {
406 
407  Qvector[js][j] = NULL;
408 
409  }
410  }
411 
412 
413 
414 } // void AliAnalysisTaskStudentsML::InitializeArrays()
415 
416 //=======================================================================================================================
417 
419 {
420  // Book and nest all lists nested in the base list fHistList.
421 
422  // a) Book and nest lists for control histograms;
423  // b) Book and nest lists for final results.
424 
425  TString sMethodName = "void AliAnalysisTaskStudentsML::BookAndNestAllLists()";
426  if(!fHistList){Fatal(sMethodName.Data(),"fHistList is NULL");}
427 
428  // a) Book and nest lists for control histograms:
430  fControlHistogramsList->SetName("ControlHistograms");
431  fControlHistogramsList->SetOwner(kTRUE);
432  fHistList->Add(fControlHistogramsList);
433 
434  // b) Book and nest lists for final results:
435  fFinalResultsList = new TList();
436  fFinalResultsList->SetName("FinalResults");
437  fFinalResultsList->SetOwner(kTRUE);
439 
440 } // void AliAnalysisTaskStudentsML::BookAndNestAllLists()
441 
442 //=======================================================================================================================
443 
445 {
446  // Book all control histograms.
447 
448  // a) Book histogram to hold pt distribution;
449  // b) Book histogram to hold phi distribution;
450  // c) Book histogram to hold eta distribution;
451  // d) Book histogram to hold multiplicty distribution;
452  // e) Book histogram to debug
453 
454 
455 
456 
457  // a) Book histogram to hold pt spectra:
458  fPtHist = new TH1F("fPtHist","atrack->Pt()",fNbins,fMinBin,fMaxBin);
459  fPtHist->SetStats(kFALSE);
460  fPtHist->SetFillColor(kBlue-10);
461  fPtHist->GetXaxis()->SetTitle("p_{t}");
463 
464  // b) Book histogram to hold phi distribution:
465  fPhiHist = new TH1F("fPhiHist","Phi Distribution",1000,0.,6.3);
466  fPhiHist->GetXaxis()->SetTitle("Phi");
467  fPhiHist->SetLineColor(4);
469 
470  // c) Book histogram to hold eta distribution:
471  fEtaHist = new TH1F("fEtaHist","Eta Distribution",1000,-1.,1.);
472  fEtaHist->GetXaxis()->SetTitle("Eta");
473  fEtaHist->SetLineColor(4);
475 
476  // d) Book histogam to hold multiplicty distribution:
477  fMultiHisto = new TH1F("fMultiHisto","Multiplicity",500,0.,500.); //histogram for multiplicity
478  fMultiHisto->GetXaxis()->SetTitle("Multiplicity M");
480 
481  // e) Book histogram to debug
482  fCounterHistogram = new TH1F("fCounterHistogram","Histogram for some checks",3,0.,3.); //histogram for multiplicity
484 
485 } //void AliAnalysisTaskStudentsML::BookControlHistograms()
486 
487 //=======================================================================================================================
488 
490 {
491  // Book all histograms to hold the final results.
492 
493  fCentrality = new TProfile("fCentrality","Result Analysis Centrality Dependence",9,0.,9.); //histogram for multiplicity
494  fCentrality->GetXaxis()->SetTitle("Centrality");
495  fCentrality->GetYaxis()->SetTitle("flow");
497  fCentrality->Sumw2();
498 
499  fCentralitySecond = new TProfile("fCentralitySecond","Result Analysis Centrality Dependence",9,0.,9.); //histogram for multiplicity
500  fCentralitySecond->GetXaxis()->SetTitle("Centrality");
501  fCentralitySecond->GetYaxis()->SetTitle("flow");
503  fCentralitySecond->Sumw2();
504 
505  Cosmetics();
506 
507  //*** for testing
508  /*func1=new TF1("func1","(1/(TMath::TwoPi()))*(1 + 2*[1]*TMath::Cos(2*(x-[2])))",0.,TMath::TwoPi()); //defining function func1 with parameter [1] for v_2 and [2] as pis_2
509  func1->SetParameter(1,0.05); //setting v_2 (aka [1]) fix to 0.05*/
510 
511 } // void AliAnalysisTaskStudentsML::BookFinalResultsHistograms()
512 
513 //=======================================================================================================================
514 
516 {
517  // Book everything here.
518 
519  for(Int_t cs=0;cs<2;cs++)
520  {
521  for(Int_t c=0;c<fMaxCorrelator;c++)
522  {
523 
524  fRecursion[cs][c] = new TProfile("","",2,0.,2.);
525  fRecursion[cs][c]->Sumw2();
526 
527  fRecursionSecond[cs][c] = new TProfile("","",2,0.,2.);
528  fRecursionSecond[cs][c]->Sumw2();
529 
530  //NOTE for fRecursion: 1.) [cs] will say if its the real (0) or imaginary part (1)
531  // 2.) [c] gives gives the kind of correlation. [n] is the (n+2)-particle correlation
532  //3.) first bin holds value of correlation (single event). Second bin holds the weight!
533 
534 
535  } // end of for(Int_t c=0;c<fMaxCorrelator;c++)
536  } // end of for(Int_t cs=0;cs<2;cs++)
537 
538 } // void Cosmetics()
539 
540 //=======================================================================================================================
541 
543 {
544  // Calculate Q-vectors.
545 
546  // a) Make sure all Q-vectors are initially zero;
547  // b) Calculate Q-vectors for available angles and weights.
548 
549  // a) Make sure all Q-vectors are initially zero:
550  for(Int_t h=0;h<kMaxHarmonic;h++)
551  {
552  for(Int_t p=0;p<kMaxPower;p++)
553  {
554  Qvector[h][p] = TComplex(0.,0.);
555  } // for(Int_t p=0;p<kMaxPower;p++)
556  } // for(Int_t h=0;h<kMaxHarmonic;h++)
557 
558  // b) Calculate Q-vectors for available angles and weights:
559  Double_t dPhi2 = 0.; // particle angle
560  Double_t wPhi = 1.; // particle weight
561  Double_t wPhiToPowerP = 1.; // particle weight raised to power p
562  for(Int_t i=0;i<fParticles;i++) // loop over particles
563  {
564  dPhi2 = fAngles->GetAt(i);
565  if(bUseWeights){wPhi = fWeights->GetAt(i);}
566  for(Int_t h=0;h<kMaxHarmonic;h++)
567  {
568  for(Int_t p=0;p<kMaxPower;p++)
569  {
570  if(bUseWeights){wPhiToPowerP = pow(wPhi,p);}
571  Qvector[h][p] += TComplex(wPhiToPowerP*TMath::Cos(h*dPhi2),wPhiToPowerP*TMath::Sin(h*dPhi2));
572  } // for(Int_t p=0;p<kMaxPower;p++)
573  } // for(Int_t h=0;h<kMaxHarmonic;h++)
574  } // for(Int_t i=0;i<fParticles;i++) // loop over particles
575 
576 } // void CalculateQvectors()
577 
578 //=======================================================================================================================
579 
580 
582 {
583  // Using the fact that Q{-n,p} = Q{n,p}^*.
584 
585  if(n>=0){return Qvector[n][p];}
586  return TComplex::Conjugate(Qvector[-n][p]);
587 
588 } // TComplex AliAnalysisTaskStudentsML::Q(Int_t n, Int_t p)
589 
590 
591 
592 
593 //========================================================================================================================
594 
595 
596 TComplex AliAnalysisTaskStudentsML::Recursion(Int_t n, Int_t* harmonic, Int_t mult = 1, Int_t skip = 0)
597 {
598  // Calculate multi-particle correlators by using recursion (an improved faster version) originally developed by
599  // Kristjan Gulbrandsen (gulbrand@nbi.dk).
600 
601  Int_t nm1 = n-1;
602  TComplex c(Q(harmonic[nm1], mult));
603  if (nm1 == 0) return c;
604  c *= Recursion(nm1, harmonic);
605  if (nm1 == skip) return c;
606 
607  Int_t multp1 = mult+1;
608  Int_t nm2 = n-2;
609  Int_t counter1 = 0;
610  Int_t hhold = harmonic[counter1];
611  harmonic[counter1] = harmonic[nm2];
612  harmonic[nm2] = hhold + harmonic[nm1];
613  TComplex c2(Recursion(nm1, harmonic, multp1, nm2));
614  Int_t counter2 = n-3;
615  while (counter2 >= skip) {
616  harmonic[nm2] = harmonic[counter1];
617  harmonic[counter1] = hhold;
618  ++counter1;
619  hhold = harmonic[counter1];
620  harmonic[counter1] = harmonic[nm2];
621  harmonic[nm2] = hhold + harmonic[nm1];
622  c2 += Recursion(nm1, harmonic, multp1, counter2);
623  --counter2;
624  }
625  harmonic[nm2] = harmonic[counter1];
626  harmonic[counter1] = hhold;
627 
628  if (mult == 1) return c-c2;
629  return c-Double_t(mult)*c2;
630 
631 } // TComplex AliFlowAnalysisWithMultiparticleCorrelations::Recursion(Int_t n, Int_t* harmonic, Int_t mult = 1, Int_t skip = 0)
632 
633 
634 
635 //========================================================================================================================
636 
637 
639 {
640 
641 
642  // Calculate Q-vectors for available angles and weights;
644 
645  // Calculate n-particle correlations from Q-vectors (using recursion):
646 
647 
648  if(kNumber==2)
649  {
650  Int_t harmonicsTwoNum[2] = {kh1,kh2};
651  Int_t harmonicsTwoDen[2] = {0,0};
652  TComplex twoRecursion = Recursion(2,harmonicsTwoNum)/Recursion(2,harmonicsTwoDen).Re();
653  Double_t wTwoRecursion = Recursion(2,harmonicsTwoDen).Re();
654  fRecursion[0][0]->Fill(0.5,twoRecursion.Re(),wTwoRecursion); // <<cos(h1*phi1+h2*phi2)>>
655  fRecursion[0][0]->Fill(1.5,wTwoRecursion,1.); //weight
656  fRecursion[1][0]->Fill(0.5,twoRecursion.Im(),wTwoRecursion); // <<sin(h1*phi1+h2*phi2)>>
657  fRecursion[1][0]->Fill(1.5,wTwoRecursion,1.); //weight
658  }// 2-p correlation
659 
660  if(kNumber==3)
661  {
662  Int_t harmonicsThreeNum[3] = {kh1,kh2,kh3};
663  Int_t harmonicsThreeDen[3] = {0,0,0};
664  TComplex threeRecursion = Recursion(3,harmonicsThreeNum)/Recursion(3,harmonicsThreeDen).Re();
665  Double_t wThreeRecursion = Recursion(3,harmonicsThreeDen).Re();
666  fRecursion[0][1]->Fill(0.5,threeRecursion.Re(),wThreeRecursion); // <<cos(h1*phi1+h2*phi2+h3*phi3)>>
667  fRecursion[0][1]->Fill(1.5,wThreeRecursion,1.); //weight
668  fRecursion[1][1]->Fill(0.5,threeRecursion.Im(),wThreeRecursion); // <<sin(h1*phi1+h2*phi2+h3*phi3)>>
669  fRecursion[1][1]->Fill(1.5,wThreeRecursion,1.); //weight
670  } // 3-p correlation
671 
672  if(kNumber==4)
673  {
674  Int_t harmonicsFourNum[4] = {kh1,kh2,kh3,kh4};
675  Int_t harmonicsFourDen[4] = {0,0,0,0};
676  TComplex fourRecursion = Recursion(4,harmonicsFourNum)/Recursion(4,harmonicsFourDen).Re();
677  Double_t wFourRecursion = Recursion(4,harmonicsFourDen).Re();
678  fRecursion[0][2]->Fill(0.5,fourRecursion.Re(),wFourRecursion); // <<cos(h1*phi1+h2*phi2+h3*phi3+h4*phi4)>>
679  fRecursion[0][2]->Fill(1.5,wFourRecursion,1.); //weigh
680  fRecursion[1][2]->Fill(0.5,fourRecursion.Im(),wFourRecursion); // <<<sin(h1*phi1+h2*phi2+h3*phi3+h4*phi4)>>
681  fRecursion[1][2]->Fill(1.5,wFourRecursion,1.); //weigh
682  }// 4-p correlation
683 
684  if(kNumber==5)
685  {
686  Int_t harmonicsFiveNum[5] = {kh1,kh2,kh3,kh4,kh5};
687  Int_t harmonicsFiveDen[5] = {0,0,0,0,0};
688  TComplex fiveRecursion = Recursion(5,harmonicsFiveNum)/Recursion(5,harmonicsFiveDen).Re();
689  Double_t wFiveRecursion = Recursion(5,harmonicsFiveDen).Re();
690  fRecursion[0][3]->Fill(0.5,fiveRecursion.Re(),wFiveRecursion); // <<cos(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5)>>
691  fRecursion[0][3]->Fill(1.5,wFiveRecursion,1.);
692  fRecursion[1][3]->Fill(0.5,fiveRecursion.Im(),wFiveRecursion); // <<<sin(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5)>>
693  fRecursion[1][3]->Fill(1.5,wFiveRecursion,1.);
694  }// 5-p correlation
695 
696  if(kNumber==6)
697  {
698  Int_t harmonicsSixNum[6] = {kh1,kh2,kh3,kh4,kh5,kh6};
699  Int_t harmonicsSixDen[6] = {0,0,0,0,0,0};
700  TComplex sixRecursion = Recursion(6,harmonicsSixNum)/Recursion(6,harmonicsSixDen).Re();
701  Double_t wSixRecursion = Recursion(6,harmonicsSixDen).Re();
702  fRecursion[0][4]->Fill(0.5,sixRecursion.Re(),wSixRecursion); // <<cos(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6)>>
703  fRecursion[0][4]->Fill(1.5,wSixRecursion,1.);
704  fRecursion[1][4]->Fill(0.5,sixRecursion.Im(),wSixRecursion); // <<<sin(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6)>>
705  fRecursion[1][4]->Fill(1.5,wSixRecursion,1.);
706 
707 
708  Int_t harmonicsSixNumSecond[6] = {ka1,ka2,ka3,ka4,ka5,ka6};
709  Int_t harmonicsSixDenSecond[6] = {0,0,0,0,0,0};
710  TComplex sixRecursionSecond = Recursion(6,harmonicsSixNumSecond)/Recursion(6,harmonicsSixDenSecond).Re();
711  Double_t wSixRecursionSecond = Recursion(6,harmonicsSixDenSecond).Re();
712 
713 
714  fRecursionSecond[0][4]->Fill(0.5,sixRecursionSecond.Re(),wSixRecursionSecond); // <<cos(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6)>>
715  fRecursionSecond[0][4]->Fill(1.5,wSixRecursionSecond,1.);
716  fRecursionSecond[1][4]->Fill(0.5,sixRecursionSecond.Im(),wSixRecursionSecond); // <<<sin(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6)>>
717  fRecursionSecond[1][4]->Fill(1.5,wSixRecursionSecond,1.);
718 
719 
720 
721  }// 6-p correlation
722 
723 
724  if(kNumber==7)
725  {
726  Int_t harmonicsSevenNum[7] = {kh1,kh2,kh3,kh4,kh5,kh6,kh7};
727  Int_t harmonicsSevenDen[7] = {0,0,0,0,0,0,0};
728  TComplex sevenRecursion = Recursion(7,harmonicsSevenNum)/Recursion(7,harmonicsSevenDen).Re();
729  Double_t wSevenRecursion = Recursion(7,harmonicsSevenDen).Re();
730  fRecursion[0][5]->Fill(0.5,sevenRecursion.Re(),wSevenRecursion); // <<cos(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6+h7*phi7)>>
731  fRecursion[0][5]->Fill(1.5,wSevenRecursion,1.);
732  fRecursion[1][5]->Fill(0.5,sevenRecursion.Im(),wSevenRecursion); // <<<sin(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6+h7*phi7)>>
733  fRecursion[1][5]->Fill(1.5,wSevenRecursion,1.);
734  }// 7-p correlation
735 
736 
737  if(kNumber==8)
738  {
739  Int_t harmonicsEightNum[8] = {kh1,kh2,kh3,kh4,kh5,kh6,kh7,kh8};
740  Int_t harmonicsEightDen[8] = {0,0,0,0,0,0,0,0};
741  TComplex eightRecursion = Recursion(8,harmonicsEightNum)/Recursion(8,harmonicsEightDen).Re();
742  Double_t wEightRecursion = Recursion(8,harmonicsEightDen).Re();
743  fRecursion[0][6]->Fill(0.5,eightRecursion.Re(),wEightRecursion); // <<cos(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6+h7*phi7+h8*phi8)>>
744  fRecursion[0][6]->Fill(1.5,wEightRecursion,1.);
745  fRecursion[1][6]->Fill(0.5,eightRecursion.Im(),wEightRecursion); // <<<sin(h1*phi1+h2*phi2+h3*phi3+h4*phi4+h5*phi5+h6*phi6+h7*phi7+h8*phi8)>>
746  fRecursion[1][6]->Fill(1.5,wEightRecursion,1.);
747  }// 8-p correlation
748 
749  if(kNumber!=2 && kNumber!=3 && kNumber!=4 && kNumber!=5 && kNumber!=6 && kNumber!=7 && kNumber!=8)
750  {
751  return;
752  }
753 
754 
755 
756  }//void Correlation()
757 
758 
759 //========================================================================================================================
760 
761 
Int_t charge
double Double_t
Definition: External.C:58
virtual void Terminate(Option_t *)
TProfile * fCentrality
//[fMaxHarmonic*fMaxCorrelator+1][fMaxCorrelator+1]
TCanvas * c
Definition: TestFitELoss.C:172
TArrayD * fWeights
Azimuthal angles.
int Int_t
Definition: External.C:63
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
TComplex Recursion(Int_t n, Int_t *harmonic, Int_t mult, Int_t skip)