Improved version, kinetree included
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowEvent.cxx
1 //////////////////////////////////////////////////////////////////////
2 //
3 // $Id$
4 //
5 // Author: Emanuele Simili
6 //
7 //////////////////////////////////////////////////////////////////////
8 //
9 //_____________________________________________________________
10 //
11 // Description: 
12 //         AliFlowEvent is the basic class to perform flow study in ALICE. 
13 // The AliFlowEvent object contains data mebers wich summarize the ESD 
14 // informations that are most useful for flow study, together with a 
15 // collection of tracks (AliFlowTrackCollection) and reconstructed v0s. 
16 // The class also implements member functions to calculate Multiplicity, 
17 // Mean P and Pt, Q (the event plane vector), Psi (the event plane angle), 
18 // normalized q (Q/sqrt(Mult) and Q/sqrt(sum of weights^2)), moreover, 
19 // functions to separate random or eta sub-events, to fill the bayesian vector 
20 // of particle abundance (for bayesian p.id calculation), to plug-in weights
21 // (to make the tracks' distribution isotropic and to enhance the resolution), 
22 // and to select the tracks that are used in the event plane calculation (for 
23 // each harmonic, selection, and subevent. 
24 // AliFlowEvent supports the standard analysis as well as the Cumulant 
25 // Analysis: method to extract the generating functions for the old and new 
26 // cumulant methods is also there.
27 // The Flow Event structure (AliFlowEvent, AliFlowTrack, AliFlowV0 and 
28 // AliFlowSelection) is independent from AliRoot, i.e. once the FlowEvents 
29 // are filled from the AliESD or the MC KineTree (see AliFlowInterface: here 
30 // AliRoot is needed), the flow analysis itself can be performed just under 
31 // ROOT (see AliFlowAnalysisMaker).
32 //
33 // AliFlowEvent is adapted from the original class StFlowEvent, succesfully 
34 // employed to study flow in the STAR experiment at RICH.
35 // Original Authors:                 Raimond Snellings & Art Poskanzer
36 //
37
38 #include "AliFlowEvent.h"
39 #include "AliFlowTrack.h"
40 #include "AliFlowV0.h"
41 #include "AliFlowSelection.h"
42 #include "AliFlowConstants.h"
43
44 #include "TVector.h"
45 #include "TVector2.h"
46 #include "TVector3.h"
47 #include "TClonesArray.h"
48 #include "TMath.h"
49 #include "TRandom.h"
50 #include "TString.h"
51 #include <iostream>
52 using namespace std; //required for resolving the 'cout' symbol
53
54 // - Flags & Sets
55 Bool_t  AliFlowEvent::fgPtWgt          = kFALSE ;  // gives pT-based weights
56 Bool_t  AliFlowEvent::fgEtaWgt         = kFALSE ;  // gives eta-based weights
57 Bool_t  AliFlowEvent::fgOnePhiWgt      = kTRUE  ;  // kTRUE --> ENABLEs SINGLE WEIGHT ISTOGRAM , kFALSE --> ENABLEs 3 WEIGHT ISTOGRAMS
58 Bool_t  AliFlowEvent::fgNoWgt          = kFALSE ;  // No Weight is used
59 // - Eta Sub-Events (later used to calculate the resolution)
60 Bool_t  AliFlowEvent::fgEtaSubs        = kFALSE ;  // makes eta subevents
61 Bool_t  AliFlowEvent::fgCustomRespFunc = kFALSE ;  // custom "detector response function" is used for P.Id (see AliFlowTrack)
62
63 ClassImp(AliFlowEvent) 
64 //-----------------------------------------------------------
65 AliFlowEvent::AliFlowEvent(Int_t lenght):
66   fTrackCollection(0x0), fV0Collection(0x0)     
67 {
68  // Default constructor: initializes the ObjArray of FlowTracks and FlowV0s, 
69  // cleans the internal variables, sets all the weights to 1, sets default flags.
70  
71  fEventID = 0 ; 
72  fRunID = 0 ;                  
73  fOrigMult = 0 ;
74  fL0Trigger = 0 ;                              
75  fCentrality = -1 ;
76  fZDCpart = 0 ;
77  for(int zz=0;zz<3;zz++) { fZDCenergy[zz] = 0. ; }
78  for(int i=0;i<AliFlowConstants::kHars;i++) { fExtPsi[i] = 0. ; fExtRes[i] = 0.; }
79  
80  // Make a new track collection
81  fTrackCollection =  new TClonesArray("AliFlowTrack",lenght) ; fTrackCollection->BypassStreamer(kTRUE) ;
82  fV0Collection    =  new TClonesArray("AliFlowV0",lenght)    ; fV0Collection->BypassStreamer(kTRUE)    ;
83  
84  // Set Weights Arrays to 1 (default)
85  for(int nS=0;nS<AliFlowConstants::kSels;nS++)
86  {
87   for(int nH=0;nH<AliFlowConstants::kHars;nH++) 
88   {
89    for(int nP=0;nP<AliFlowConstants::kPhiBins;nP++) 
90    { 
91     // enable this with: SetOnePhiWgt()  
92     fPhiWgt[nS][nH][nP] = 1.      ; // cout << nS << nH << nP << "  val  " << fPhiWgt[nS][nH][nP] << endl ; 
93     // enable these with: SetFirstLastPhiWgt()  
94     fPhiWgtPlus[nS][nH][nP]  = 1. ; // cout << nS << nH << nP << "  val  " << fPhiWgtPlus[nS][nH][nP] << endl ; 
95     fPhiWgtMinus[nS][nH][nP] = 1. ; // cout << nS << nH << nP << "  val  " << fPhiWgtMinus[nS][nH][nP] << endl ; 
96     fPhiWgtCross[nS][nH][nP] = 1. ; // cout << nS << nH << nP << "  val  " << fPhiWgtCross[nS][nH][nP] << endl ; 
97    }
98   }
99  }
100  //for(int nH=0;nH<AliFlowConstants::kHars;nH++) { fExtPsi[nH] = 0. ; fExtRes[nH] = 0. ; }
101  
102  // The Expected particles abundance is taken directly from AliFlowConstants::fgBayesian[] (see Bayesian P.Id.)
103  
104  fDone = kFALSE ;
105 }
106 //-----------------------------------------------------------
107 AliFlowEvent::~AliFlowEvent() 
108 {
109  // Default distructor: deletes the ObjArrays. 
110  
111  fTrackCollection->Delete() ; delete fTrackCollection ;
112  fV0Collection->Delete()    ; delete fV0Collection ;
113 }
114 //-------------------------------------------------------------
115
116 //-------------------------------------------------------------
117 Double_t AliFlowEvent::PhiWeightRaw(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const  
118 {
119  // Weight for making the event plane isotropic in the lab.
120
121  float phi = pFlowTrack->Phi() ; if(phi < 0.) { phi += 2*TMath::Pi() ; }
122  Double_t  eta = (Double_t)pFlowTrack->Eta() ;
123  int n = (int)((phi/(2*TMath::Pi()))*AliFlowConstants::kPhiBins);
124
125  Double_t phiWgt = 1. ;
126  if(OnePhiWgt()) 
127  {
128   phiWgt = (Double_t)fPhiWgt[selN][harN][n]; //cout << "Selection " << selN << " ; Harmonic " << harN << " ; PhiBin " << n << "  -  Wgt = " << phiWgt << endl ; 
129  } 
130  else if(FirstLastPhiWgt())
131  {
132   float zFirstPoint = pFlowTrack->ZFirstPoint(); // float zLastPoint = pFlowTrack->ZLastPoint();
133   
134   if (zFirstPoint > 0. && eta > 0.)     { phiWgt = (Double_t)fPhiWgtPlus[selN][harN][n] ; } 
135   else if(zFirstPoint < 0. && eta < 0.) { phiWgt = (Double_t)fPhiWgtMinus[selN][harN][n] ; } 
136   else                                  { phiWgt = (Double_t)fPhiWgtCross[selN][harN][n] ; }
137  } 
138
139  return phiWgt ;
140 }
141 //-------------------------------------------------------------
142 Double_t AliFlowEvent::Weight(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const  
143
144  // Weight for enhancing the resolution (eta gives sign +/- for Odd Harmonics)
145
146  if(selN>AliFlowConstants::kSels) { selN = 0 ; }
147  bool oddHar = (harN+1) % 2 ;
148  Double_t phiWgt = 1. ;
149  if(PtWgt()) 
150  {
151   Double_t pt = (Double_t)pFlowTrack->Pt();
152   if(pt < AliFlowConstants::fgPtWgtSaturation) { phiWgt *= pt ; } 
153   else                            { phiWgt *= AliFlowConstants::fgPtWgtSaturation ; } // pt weighting going constant
154  }
155  Double_t eta = (Double_t)pFlowTrack->Eta();
156  Double_t etaAbs = TMath::Abs(eta);
157  if(EtaWgt() && oddHar)     { phiWgt *= etaAbs ; }
158  if(oddHar && eta < 0.)     { phiWgt *= -1. ; }
159
160  return phiWgt ;
161 }
162 //-------------------------------------------------------------
163 Double_t AliFlowEvent::PhiWeight(Int_t selN, Int_t harN, AliFlowTrack* pFlowTrack) const
164 {
165  // Weight for making the event plane isotropic in the lab and enhancing the resolution
166  // (it simply rerurns PhiWeightRaw() * Weight()). If fgNoWgt = kTRUE, returns +/-1 ,
167  // basing on Sign(eta), for odd harmonics .
168  
169  if(fgNoWgt) // no weights (but +/- according to eta)
170  { 
171   bool oddHar = (harN+1) % 2 ;
172   if(oddHar) { return TMath::Sign((Double_t)1.,(Double_t)pFlowTrack->Eta()) ; }
173   else       { return (Double_t)1. ; }
174  } 
175  Double_t phiWgtRaw = PhiWeightRaw(selN, harN, pFlowTrack);
176  Double_t weight = Weight(selN, harN, pFlowTrack);
177  if(AliFlowConstants::fgDebug) { cout << "[PhiWeight]: phiWgtRaw = " << phiWgtRaw << " , weight = " << weight << " , eta = " << pFlowTrack->Eta() << endl ; }
178
179  return phiWgtRaw * weight;
180 }
181 //-------------------------------------------------------------
182 Int_t AliFlowEvent::Mult(AliFlowSelection* pFlowSelect) const   
183 {
184  // Multiplicity of tracks in the specified Selection 
185  
186  if(fDone) 
187  {
188   int sub = pFlowSelect->Sub() ;
189   if(sub<0) { return fMult[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
190   else      { return fMultSub[sub][pFlowSelect->Sel()][pFlowSelect->Har()] ; }
191  }
192  // -
193  Int_t mult = 0;
194  Int_t itr ;
195  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
196  {
197   AliFlowTrack* pFlowTrack ;
198   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
199   if(pFlowSelect->Select(pFlowTrack)) { mult++ ; }
200  }
201  return mult;
202 }
203 //-------------------------------------------------------------
204 Float_t AliFlowEvent::MeanPt(AliFlowSelection* pFlowSelect) const
205 {
206  // Mean pt of tracks in the specified Selection 
207
208  Double_t meanPt = 0. ;
209  Float_t sumPt = 0. ;
210  UInt_t mult = 0 ;
211  Int_t itr ;
212  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
213  {
214   AliFlowTrack* pFlowTrack ;
215   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
216   if (pFlowSelect->Select(pFlowTrack)) 
217   {
218    sumPt += pFlowTrack->Pt();
219    mult++;
220   }
221  }
222  if(mult) { meanPt = sumPt/(float)mult ; }
223  
224  return meanPt ;
225 }
226 //-------------------------------------------------------------
227 TVector2 AliFlowEvent::Q(AliFlowSelection* pFlowSelect) const
228 {
229  // Event plane vector for the specified Selection 
230
231  if(fDone) 
232  { 
233   int sub = pFlowSelect->Sub() ;
234   if(sub<0) { return fQ[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
235   else      { return fQSub[sub][pFlowSelect->Sel()][pFlowSelect->Har()] ; }
236  } 
237  // -
238  TVector2 mQ ;
239  Double_t mQx=0. , mQy=0. ;
240  int    selN  = pFlowSelect->Sel() ;
241  int    harN  = pFlowSelect->Har() ;
242  double order = (double)(harN + 1) ;
243
244  Int_t itr ;
245  for(itr=0;itr<TrackCollection()->GetEntries();itr++)
246  {
247   AliFlowTrack* pFlowTrack ;
248   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
249   if(pFlowSelect->Select(pFlowTrack)) 
250   {
251    double phiWgt = PhiWeight(selN, harN, pFlowTrack);  
252    float phi = pFlowTrack->Phi();                      
253    mQx += phiWgt * cos(phi * order) ;
254    mQy += phiWgt * sin(phi * order) ;
255    if(AliFlowConstants::fgDebug) { cout << itr << " phi = " << phi << " ,  wgt = " << phiWgt << endl ; }
256   }
257  }
258  mQ.Set(mQx, mQy);
259
260  return mQ;
261 }
262 //-------------------------------------------------------------
263 Float_t AliFlowEvent::Psi(AliFlowSelection* pFlowSelect) const
264 {
265  // Event plane angle for the specified Selection 
266
267  int   harN = pFlowSelect->Har() ;
268  float order = (float)(harN + 1) ;
269  Float_t psi = 0. ;
270
271  TVector2 mQ = Q(pFlowSelect);
272  if(mQ.Mod())  // if vector is not 0
273  {
274   psi= mQ.Phi() / order ;
275   if (psi < 0.) { psi += 2*TMath::Pi() / order ; }
276  }
277  return psi;
278 }
279 //-------------------------------------------------------------
280 TVector2 AliFlowEvent::NormQ(AliFlowSelection* pFlowSelect) const       
281 {
282  // Normalized Q = Q/sqrt(sum of weights^2) for the specified Selection 
283
284  if(fDone) 
285  {
286   TVector2 mQ = fQ[pFlowSelect->Sel()][pFlowSelect->Har()] ;   
287   double sumOfWeightSqr = fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ;
288   if(sumOfWeightSqr) { mQ /= TMath::Sqrt(sumOfWeightSqr) ; }
289   else { mQ.Set(0.,0.) ; }
290   return mQ ;
291  }
292  // -
293  TVector2 mQ ;
294  Double_t mQx=0. , mQy=0. ;
295  double sumOfWeightSqr = 0 ;
296  int selN     = pFlowSelect->Sel() ;
297  int harN     = pFlowSelect->Har() ;
298  double order = (double)(harN + 1) ;
299  Int_t itr ;
300  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
301  {
302   AliFlowTrack* pFlowTrack ;
303   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
304   if (pFlowSelect->Select(pFlowTrack)) 
305   {
306    double phiWgt = PhiWeight(selN, harN, pFlowTrack);
307    sumOfWeightSqr += phiWgt*phiWgt;
308
309    float phi = pFlowTrack->Phi();
310    mQx += phiWgt * cos(phi * order);
311    mQy += phiWgt * sin(phi * order);
312   }
313  }
314  if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
315  else { mQ.Set(0.,0.); }
316   
317  return mQ;
318 }
319 //-------------------------------------------------------------
320 Float_t AliFlowEvent::OldQ(AliFlowSelection* pFlowSelect) const 
321
322  // Magnitude of normalized Q vector (without pt or eta weighting) for the specified Selection 
323
324  TVector2 mQ ;
325  Double_t mQx = 0. , mQy = 0. ;
326  int selN     = pFlowSelect->Sel() ;
327  int harN     = pFlowSelect->Har() ;
328  double order = (double)(harN + 1) ;
329  double sumOfWeightSqr = 0 ;
330
331  Int_t itr ;
332  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
333  {
334   AliFlowTrack* pFlowTrack ;
335   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
336   if(pFlowSelect->Select(pFlowTrack)) 
337   {
338    double phiWgt = PhiWeightRaw(selN, harN, pFlowTrack); // Raw
339    sumOfWeightSqr += phiWgt*phiWgt;
340    float phi = pFlowTrack->Phi();
341    mQx += phiWgt * cos(phi * order);
342    mQy += phiWgt * sin(phi * order);
343   }
344  }
345  if(sumOfWeightSqr) { mQ.Set(mQx/TMath::Sqrt(sumOfWeightSqr), mQy/TMath::Sqrt(sumOfWeightSqr)); }
346  else { mQ.Set(0.,0.); }
347   
348  return mQ.Mod();
349 }
350 //-----------------------------------------------------------------------
351 Double_t AliFlowEvent::NewG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
352
353  // Generating function for the new cumulant method. Eq. 3 in the Practical Guide 
354
355  int selN     = pFlowSelect->Sel();
356  int harN     = pFlowSelect->Har();
357  double order = (double)(harN + 1);
358
359  double mult = (double)Mult(pFlowSelect);
360  Double_t theG = 1.;
361
362  Int_t itr ;
363  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
364  {
365   AliFlowTrack* pFlowTrack ;
366   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
367   if (pFlowSelect->Select(pFlowTrack)) 
368   {
369    double phiWgt = PhiWeight(selN, harN, pFlowTrack);      
370    float phi = pFlowTrack->Phi();
371    theG *= (1. + (phiWgt/mult) * (2.* Zx * cos(phi * order) + 2.* Zy * sin(phi * order) ) );            
372   }
373  }
374  return theG;
375 }
376 //-----------------------------------------------------------------------
377 Double_t AliFlowEvent::OldG(AliFlowSelection* pFlowSelect, Double_t Zx, Double_t Zy) const
378
379  // Generating function for the old cumulant method (if expanded in Taylor
380  // series, one recovers NewG() in new new cumulant method)
381
382  TVector2 normQ = NormQ(pFlowSelect) ;
383
384  return exp(2*Zx*normQ.X() + 2*Zy*normQ.Y());
385 }
386 //-----------------------------------------------------------------------
387 Double_t AliFlowEvent::SumWeightSquare(AliFlowSelection* pFlowSelect) const
388 {
389  // Return sum of weights^2 for the specified Selection (used for normalization)
390
391  if(fDone) { return fSumOfWeightSqr[pFlowSelect->Sel()][pFlowSelect->Har()] ; }
392
393  int selN = pFlowSelect->Sel();
394  int harN = pFlowSelect->Har();
395  Double_t sumOfWeightSqr = 0;
396  Int_t itr ;
397  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
398  {
399   AliFlowTrack* pFlowTrack ;
400   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
401   if (pFlowSelect->Select(pFlowTrack)) 
402   {
403    double phiWgt = PhiWeight(selN, harN, pFlowTrack);
404    sumOfWeightSqr += phiWgt*phiWgt;
405   }
406  }
407  if(sumOfWeightSqr < 0.) { return Mult(pFlowSelect) ; }
408
409  return sumOfWeightSqr;
410 }
411 //-------------------------------------------------------------
412 Double_t AliFlowEvent::WgtMultQ4(AliFlowSelection* pFlowSelect) const 
413
414  // Used only for the old cumulant method, for getting q4 when weight is on.
415  // Replace multiplicity in Eq.(74b) by this quantity when weight is on.
416  // This is derived based on (A4) in the old cumulant paper.
417
418  int selN = pFlowSelect->Sel();
419  int harN = pFlowSelect->Har();
420  double theMult        = 0.;
421  double theMeanWj4     = 0.;
422  double theMeanWj2     = 0.;
423  double theSumOfWgtSqr = 0;
424  double phiWgtSq;
425
426  Int_t itr ;
427  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
428  {
429   AliFlowTrack* pFlowTrack ;
430   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
431   if (pFlowSelect->Select(pFlowTrack)) 
432   {
433    double phiWgt   = PhiWeight(selN, harN, pFlowTrack);
434    phiWgtSq        = phiWgt*phiWgt;
435    theSumOfWgtSqr += phiWgtSq;
436    theMeanWj4     += phiWgtSq*phiWgtSq;
437    theMult        += 1.;      
438   }
439  }
440  if (theMult <= 0.) return theMult;
441
442  theMeanWj4 /= theMult;
443  theMeanWj2  = theSumOfWgtSqr / theMult;
444
445  return (theSumOfWgtSqr*theSumOfWgtSqr)/(theMult*(-theMeanWj4+2*theMeanWj2*theMeanWj2));
446 }
447 //-------------------------------------------------------------
448 Double_t AliFlowEvent::WgtMultQ6(AliFlowSelection* pFlowSelect) const 
449
450  // Used only for the old cumulant method. For getting q6 when weight is on.
451  // Replace multiplicity in Eq.(74c) by this quantity when weight is on.
452  // This is derived based on (A4) in the old cumulant paper.
453
454  int selN = pFlowSelect->Sel();
455  int harN = pFlowSelect->Har();
456  double theMult        = 0.;
457  double theMeanWj6     = 0.;
458  double theMeanWj4     = 0.;
459  double theMeanWj2     = 0.;
460  double theSumOfWgtSqr = 0;
461  double phiWgtSq;
462
463  Int_t itr ;
464  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
465  {
466   AliFlowTrack* pFlowTrack ;
467   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
468   if (pFlowSelect->Select(pFlowTrack)) 
469   {
470    double phiWgt   = PhiWeight(selN, harN, pFlowTrack);
471    phiWgtSq        = phiWgt*phiWgt;
472    theSumOfWgtSqr += phiWgtSq;
473    theMeanWj4     += phiWgtSq*phiWgtSq;
474    theMeanWj6     += phiWgtSq*phiWgtSq*phiWgtSq;
475    theMult        += 1.;
476   }
477  } 
478  if (theMult <= 0.) return theMult*theMult;
479
480  theMeanWj6 /= theMult;
481  theMeanWj4 /= theMult;
482  theMeanWj2  = theSumOfWgtSqr / theMult;
483
484  return 4.*(theSumOfWgtSqr*theSumOfWgtSqr*theSumOfWgtSqr)/(theMult*(theMeanWj6-9.*theMeanWj2*theMeanWj4+12.*theMeanWj2*theMeanWj2*theMeanWj2));
485 }
486 //-------------------------------------------------------------
487 void AliFlowEvent::SetSelections(AliFlowSelection* pFlowSelect)          
488 {
489  // Sets the selection of tracks used in the Reaction Plane calculation 
490  // for the specific Harmonic and Selection - this does not cut trow away 
491  // tracks from the event, neither exclude them from flow study. See also 
492  // the AliFlowSelection class.
493  // Strategy of Selection: 
494  //                     For the specific Harmonic and Selection, IF cuts 
495  // are defined (such that low<high) and IF the track satisfies them, THEN 
496  // the respective track's flag (bool AliFlowTrack::mSelection[har][sel]) 
497  // is set kTRUE so that the track, from now on, will be included in the 
498  // R.P. determination for that selection.
499  // If NO cuts are defined -> ALL Flags are setted kTRUE (all particle 
500  // used for all the Reaction Planes -> no difference in Psi[har][sel]).
501  // -------------------------------------------------------------------
502  // The first selection (all harmonics) is set kTRUE : no conditions.
503
504  Int_t itr ;
505  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
506  {
507   AliFlowTrack* pFlowTrack ;
508   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
509   pFlowTrack->ResetSelection() ;                // this re-sets all the mSelection flags to 0
510
511  // * this sets all the selection n.[0] flag kTRUE (all harmonics) *
512   for(int harN=0;harN<AliFlowConstants::kHars;harN++) { pFlowTrack->SetSelect(harN,0) ; } 
513
514  // Track need to be Constrainable
515   if(pFlowSelect->ConstrainCut() && !pFlowTrack->IsConstrainable()) continue ;
516
517  // PID - gets rid of the track with AliFlowTrack::Pid() != AliFlowSelection::Pid() (if there)
518   if(pFlowSelect->Pid()[0] != '\0') 
519   {
520    if(strstr(pFlowSelect->Pid(), "h")!=0) 
521    {
522     int charge = pFlowTrack->Charge() ;
523     if(strcmp("h+",pFlowSelect->Pid())==0 && charge != 1)  continue;
524     if(strcmp("h-",pFlowSelect->Pid())==0 && charge != -1) continue;
525    } 
526    else 
527    {
528     Char_t pid[10];
529     strcpy(pid, pFlowTrack->Pid());
530     if(strstr(pid, pFlowSelect->Pid())==0) continue;
531    }
532   }
533   double eta = (double)(pFlowTrack->Eta());      
534   float  pt  = pFlowTrack->Pt();                 
535   float  gDca = pFlowTrack->Dca() ;                     
536
537  // Global DCA - gets rid of the track with DCA outside the range
538   if((pFlowSelect->DcaGlobalCutHi()>pFlowSelect->DcaGlobalCutLo()) && (gDca<pFlowSelect->DcaGlobalCutLo() || gDca>pFlowSelect->DcaGlobalCutHi())) continue ;
539
540  // Pt & Eta - this is done differently for different Harmonic & Selection
541   for(int selN = 1; selN < AliFlowConstants::kSels; selN++)               // not even consider the 0th selection (no cut applied there)
542   {
543   // min. TPC hits required
544    if(pFlowSelect->NhitsCut(selN) && (pFlowTrack->FitPtsTPC()<pFlowSelect->NhitsCut(selN))) continue ;
545
546    for(int harN = 0; harN < AliFlowConstants::kHars; harN++) 
547    {
548    // Eta - gets rid of the track with Eta outside the range
549     if((pFlowSelect->EtaCutHi(harN%2,selN)>pFlowSelect->EtaCutLo(harN%2,selN)) && (TMath::Abs(eta)<pFlowSelect->EtaCutLo(harN%2,selN) || TMath::Abs(eta)>pFlowSelect->EtaCutHi(harN%2,selN))) continue ; 
550    // Pt - gets rid of the track with Pt outside the range
551     if((pFlowSelect->PtCutHi(harN%2,selN)>pFlowSelect->PtCutLo(harN%2,selN)) && (pt<pFlowSelect->PtCutLo(harN%2,selN) || pt>pFlowSelect->PtCutHi(harN%2,selN))) continue ; 
552   
553     pFlowTrack->SetSelect(harN, selN) ;  // if cuts defined (low<high) && track is in the range -> Set [har][sel] Flag ON
554
555     if(AliFlowConstants::fgDebug) 
556     {
557      if(harN==1)
558      {
559       cout << " n. " << itr << " , pFlowTrack->PrintSelection() = " ; pFlowTrack->PrintSelection() ;
560       if(pFlowSelect->Pid()[0] != '\0') { cout << " track:  pid  " << pFlowTrack->Pid() << " = "<< pFlowSelect->Pid() << endl ; } 
561       cout << " track:  dca  " << pFlowSelect->DcaGlobalCutLo() << " < " << gDca << " < " << pFlowSelect->DcaGlobalCutHi() << endl ;
562       cout << " track:  eta  " << pFlowSelect->EtaCutLo(harN,selN) << " < |" << eta << "| < " << pFlowSelect->EtaCutHi(harN,selN) << endl ;
563       cout << " track:  pT   " << pFlowSelect->PtCutLo(harN,selN) << " < " << pt << " < " << pFlowSelect->PtCutHi(harN,selN) << endl ;
564       cout << " pFlowTrack->PrintSelection() = " ; pFlowTrack->PrintSelection() ;
565      }
566      // cout << " selN " << selN << " ,  harN " << harN%2 << " - si" << endl ;
567     }
568    }
569   }
570  }
571 }
572 //-------------------------------------------------------------
573 void AliFlowEvent::SetPids()
574 {
575  // Re-sets the tracks P.id. (using the current AliFlowConstants::fgBayesian[] array).
576  // To use the Raw P.Id (just the detector response function), set fgBayesian[] = {1,1,1,1,1,1}.
577  
578  const Int_t kCode[] = {11,13,211,321,2212,10010020} ;
579  for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++) 
580  {
581   AliFlowTrack* pFlowTrack ;
582   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
583
584   Float_t bayPid[AliFlowConstants::kPid] ;
585   if(fgCustomRespFunc)  { pFlowTrack->PidProbsC(bayPid) ; }
586   else                  { pFlowTrack->PidProbs(bayPid)  ; }
587
588   Int_t maxN = 2 ;                // if No id. -> then is a Pi
589   Float_t pidMax = bayPid[2] ;    // (if all equal, Pi probability was the first)
590   for(Int_t nP=0;nP<AliFlowConstants::kPid;nP++)
591   {
592    if(bayPid[nP]>pidMax) { maxN = nP ; pidMax = bayPid[nP] ; }
593   }
594   Int_t pdgCode = TMath::Sign(kCode[maxN],pFlowTrack->Charge()) ;     
595   pFlowTrack->SetMostLikelihoodPID(pdgCode) ;                   
596  }
597 }
598 //-------------------------------------------------------------
599 void AliFlowEvent::MakeSubEvents() const
600 {
601  // Make random or eta sub-events
602  
603  if(fgEtaSubs)  { MakeEtaSubEvents() ; }
604  else           { MakeRndSubEvents() ; }
605
606 //-------------------------------------------------------------
607 void AliFlowEvent::MakeRndSubEvents() const 
608 {
609  // Make random subevents
610  
611  int eventMult[AliFlowConstants::kHars][AliFlowConstants::kSels] = {{0}};
612  int harN, selN, subN = 0;
613  
614  // loop to count the total number of tracks for each selection
615  for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++) 
616  {
617   AliFlowTrack* pFlowTrack ;
618   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
619   for (selN = 0; selN < AliFlowConstants::kSels; selN++) 
620   {
621    for (harN = 0; harN < AliFlowConstants::kHars; harN++) 
622    {
623     if(pFlowTrack->Select(harN, selN)) { eventMult[harN][selN]++ ; }
624    }
625   }
626  }
627  // loop to set the SubEvent member
628  for (selN = 0; selN < AliFlowConstants::kSels; selN++) 
629  {
630   for (harN = 0; harN < AliFlowConstants::kHars; harN++) 
631   {
632    int subEventMult = eventMult[harN][selN] / AliFlowConstants::kSubs;
633    if (subEventMult) 
634    {
635     subN = 0;
636     int countN = 0;
637     for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++) 
638     {
639      AliFlowTrack* pFlowTrack ;
640      pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
641      if(pFlowTrack->Select(harN, selN)) 
642      {
643       pFlowTrack->SetSubevent(harN, selN, subN);
644       countN++;
645       if (countN % subEventMult == 0.) { subN++ ; }
646      }
647     }
648    }
649   }
650  }
651  return ; 
652 }
653 //-------------------------------------------------------------
654 void AliFlowEvent::MakeEtaSubEvents() const
655 {
656  // Make subevents for positive and negative eta 
657  // (when done, fgEtaSubs flag setted to kTRUE).
658  
659  int harN, selN = 0;
660  // loop to set the SubEvent member
661  for (selN = 0; selN < AliFlowConstants::kSels; selN++) 
662  {
663   for (harN = 0; harN < AliFlowConstants::kHars; harN++) 
664   {
665    for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++) 
666    {
667     AliFlowTrack* pFlowTrack ;
668     pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
669     if(pFlowTrack->Select(harN, selN)) 
670     {
671      float eta = pFlowTrack->Eta();
672      if (eta > 0.) { pFlowTrack->SetSubevent(harN, selN, 0) ; } 
673      else          { pFlowTrack->SetSubevent(harN, selN, 1) ; }
674     }
675    }
676   }
677  }
678 }
679 //-------------------------------------------------------------
680 void AliFlowEvent::RandomShuffle() 
681 {
682  // Randomly re-shuffles the tracks in the array; if a track is not
683  // primary, the reference carried by the reconstructed mother (in 
684  // the v0 array) is changed accordingly.
685
686  return ; // ...at the moment is disabled ... // 
687  
688 //  Int_t tot = 0 ;
689 //  UInt_t imax = fTrackCollection->GetEntries() ;
690 //  TRandom* rnd = new TRandom(0) ; // SetSeed(0) ;
691 //  TClonesArray* newTrackCollection = new TClonesArray("AliFlowTrack",imax) ;
692 // 
693 //  // re-arranges the ObjArray (TrackCollection())
694 //  for(UInt_t itr=0;itr<imax;itr++)
695 //  {
696 //   AliFlowTrack* pFlowTrack ;
697 //   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
698 //   
699 //   UInt_t rndNumber = rnd->Integer(imax) ;
700 //   Bool_t put = kFALSE ;
701 //   while(!put)
702 //   { 
703 //    if(!newTrackCollection->AddrAt(rndNumber)) 
704 //    { 
705 //     ... new(newTrackCollection[rndNumber]) AliFlowTrack(*pFlowTrack)  ; 
706 //     put = kTRUE ; tot++ ;            
707 //     if(AliFlowConstants::fgDebug) { cout << "  " << itr << " --> " << rndNumber << endl ; } 
708 //    }
709 //    else 
710 //    {
711 //     rndNumber++ ; if(rndNumber>=imax) { rndNumber -= imax ; }
712 //    }
713 //   }
714 //  }
715 //  if(AliFlowConstants::fgDebug) { cout << "* RandomShuffle() :  " << tot << "/" << imax << " flow tracks have been shuffled " << endl ; }  
716 //  fTrackCollection = newTrackCollection ;
717 }
718 //-----------------------------------------------------------------------
719 UInt_t AliFlowEvent::Centrality() 
720 {
721  // Returns the Centrality Class as stored
722
723  if(fCentrality < 0)  { SetCentrality() ; } 
724  return fCentrality ;
725 }
726 //-----------------------------------------------------------------------
727 void AliFlowEvent::SetCentrality() 
728 {
729  // Sets the Centrality Classes basing on Multiplicity at mid rapidity, 
730  // ... (ZDC information can be added) .
731  
732  Float_t* cent ;
733  Int_t  tracks = MultEta() ;
734
735  if(RunID() == -1)
736  { 
737   cent = AliFlowConstants::fgCentNorm ;
738   //if centrality classes are not defined, does it now (with CentNorm & MaxMult)
739   if(cent[AliFlowConstants::kCents-1] <= 1) 
740   {
741    for(Int_t ic=0;ic<AliFlowConstants::kCents;ic++)
742    {
743     cent[ic] *= AliFlowConstants::fgMaxMult ; 
744     if(AliFlowConstants::fgDebug) { cout << "Centrality[" << ic << "] = " << cent[ic] << " . " << endl ; }
745    }
746   }
747  }
748  else if((RunID() != -1) && (CenterOfMassEnergy() == 5500.))
749  {
750   cent = (Float_t*)AliFlowConstants::fgCent0 ;
751  } 
752  else // other definition of centrality are possible ...
753  {
754   cent = (Float_t*)AliFlowConstants::fgCent0 ;
755  } 
756  if      (tracks < cent[0])  { fCentrality = 0; }
757  else if (tracks < cent[1])  { fCentrality = 1; }
758  else if (tracks < cent[2])  { fCentrality = 2; }
759  else if (tracks < cent[3])  { fCentrality = 3; }
760  else if (tracks < cent[4])  { fCentrality = 4; }
761  else if (tracks < cent[5])  { fCentrality = 5; }
762  else if (tracks < cent[6])  { fCentrality = 6; }
763  else if (tracks < cent[7])  { fCentrality = 7; }
764  else if (tracks < cent[8])  { fCentrality = 8; }
765  else                        { fCentrality = 9; }
766
767  if(AliFlowConstants::fgDebug) { cout << " * Centrality Class :  " << fCentrality << " . " << endl ; }
768 }
769 //-----------------------------------------------------------------------
770 TVector AliFlowEvent::Bayesian() const 
771
772  // Returns bayesian array of particle abundances (from AliFlowConstants::)
773  
774  TVector bayes(AliFlowConstants::kPid) ; 
775  for(Int_t i=0;i<AliFlowConstants::kPid;i++) { bayes[i] = AliFlowConstants::fgBayesian[i] ; }
776  return bayes ;
777 }
778 //-----------------------------------------------------------------------
779 void AliFlowEvent::PrintFlagList() const 
780 {
781  // Prints the list of selection cuts ( called in AliFlowInterface::Finish() )
782  
783  cout << "#######################################################" << endl;
784  cout << "# Weighting and Striping:" << endl;
785  if(PtWgt()) 
786  {
787   cout << "#    PtWgt =  kTRUE " << endl ;      // (also for output of PhiWgt file?)
788   cout << "#    PtWgt Saturation =  " << AliFlowConstants::fgPtWgtSaturation << endl;
789  } 
790  else 
791  {
792   cout << "#    PtWgt = kFALSE" << endl;
793  }
794  if(EtaWgt()) 
795  {
796   cout << "#    EtaWgt = kTRUE " << endl ;      // (also for output of PhiWgt file for odd harmonics?)
797  } 
798  else 
799  {
800   cout << "#    EtaWgt = kFALSE" << endl;
801  }
802  cout << "#######################################################" << endl;  
803 }
804 //-----------------------------------------------------------------------
805 Int_t AliFlowEvent::MultEta() const
806 {
807  // Returns the multiplicity in the interval |eta|<(AliFlowConstants::fgEetaMid), used 
808  // for centrality measurement (see centrality classes in fCentrality) .
809  
810  Int_t goodtracks = 0 ;
811  for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++) 
812  {
813   AliFlowTrack* pFlowTrack ;
814   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
815   if((pFlowTrack->Charge()) && (TMath::Abs(pFlowTrack->Eta())<AliFlowConstants::fgEtaMid)) { goodtracks++ ; }
816  }
817  return goodtracks ; 
818 }
819 //-----------------------------------------------------------------------
820 Int_t AliFlowEvent::UncorrNegMult(Float_t eta)  const
821
822  // Negative multiplicity in the interval (-eta..eta)
823  // (default is  AliFlowConstants::fgEetaGood = 0.9)
824  
825  Int_t negMult = 0 ;
826  for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++) 
827  {
828   AliFlowTrack* pFlowTrack ;
829   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
830   if((pFlowTrack->Charge()<0) && (TMath::Abs(pFlowTrack->Eta())<TMath::Abs(eta))) { negMult++ ; }
831   delete pFlowTrack ;
832  }
833  return negMult; 
834 }
835 //-----------------------------------------------------------------------
836 Int_t AliFlowEvent::UncorrPosMult(Float_t eta)  const
837
838  // Positive multiplicity in the interval (-eta..eta)
839  // (default is  AliFlowConstants::fgEetaGood = 0.9)
840  
841  Int_t posMult = 0 ;
842  for(Int_t itr=0;itr<TrackCollection()->GetEntries();itr++) 
843  {
844   AliFlowTrack* pFlowTrack ;
845   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
846   if((pFlowTrack->Charge()>0) && (TMath::Abs(pFlowTrack->Eta())<TMath::Abs(eta))) { posMult++ ; }
847   delete pFlowTrack ;
848  }
849  return posMult; 
850 }
851 //-----------------------------------------------------------------------
852 TVector3 AliFlowEvent::VertexPos() const
853 {
854  // Returns primary vertex position as a TVector3
855
856  Float_t vertex[3] ;
857  VertexPos(vertex) ;
858  return TVector3(vertex) ;
859 }
860 //-----------------------------------------------------------------------
861 void AliFlowEvent::MakeAll()
862
863  // calculates all quantities in 1 shoot 
864
865  Double_t mQx[AliFlowConstants::kSels][AliFlowConstants::kHars] ;
866  Double_t mQy[AliFlowConstants::kSels][AliFlowConstants::kHars] ;
867  Double_t mQxSub[AliFlowConstants::kSubs][AliFlowConstants::kSels][AliFlowConstants::kHars] ;
868  Double_t mQySub[AliFlowConstants::kSubs][AliFlowConstants::kSels][AliFlowConstants::kHars] ;
869  // -
870  int selN, harN, subN ;
871  for(selN=0;selN<AliFlowConstants::kSels;selN++) 
872  {
873   for(harN=0;harN<AliFlowConstants::kHars;harN++) 
874   {
875    mQx[selN][harN]    = 0. ;    
876    mQy[selN][harN]    = 0. ;    
877    fMult[selN][harN]  = 0 ;     
878    fSumOfWeightSqr[selN][harN] = 0. ;
879    for(subN=0;subN<AliFlowConstants::kSubs;subN++)
880    {
881     mQxSub[subN][selN][harN]   = 0. ;
882     mQySub[subN][selN][harN]   = 0. ;
883     fMultSub[subN][selN][harN] = 0 ;
884    }
885   }
886  }
887  
888  double order = 0.  ;
889  double phiWgt = 0. ;  
890  float  phi = 0.    ;                
891  // -
892  int itr ;
893  for(itr=0;itr<TrackCollection()->GetEntries();itr++) 
894  {
895   AliFlowTrack* pFlowTrack ;
896   pFlowTrack = (AliFlowTrack*)TrackCollection()->At(itr) ;
897   phi = pFlowTrack->Phi();
898   for(selN=0;selN<AliFlowConstants::kSels;selN++) 
899   {
900    for(harN=0;harN<AliFlowConstants::kHars;harN++) 
901    {
902     order = (double)(harN+1) ;
903     if(pFlowTrack->Select(harN,selN)) 
904     {
905      phiWgt = PhiWeight(selN,harN,pFlowTrack) ;  
906      fSumOfWeightSqr[selN][harN] += phiWgt*phiWgt ;
907      mQx[selN][harN] += phiWgt * cos(phi * order) ;
908      mQy[selN][harN] += phiWgt * sin(phi * order) ;
909      fMult[selN][harN]++ ;      
910      for(subN=0;subN<AliFlowConstants::kSubs;subN++)
911      {
912       if(pFlowTrack->Select(harN,selN,subN))
913       {
914        mQxSub[subN][selN][harN] += phiWgt * cos(phi * order) ;
915        mQySub[subN][selN][harN] += phiWgt * sin(phi * order) ;
916        fMultSub[subN][selN][harN]++  ;
917       }
918      }  // sub
919     }  // if
920    }  // har
921   }  // sel
922  }  //itr
923
924  for(selN=0;selN<AliFlowConstants::kSels;selN++)  
925  {
926   for(harN=0;harN<AliFlowConstants::kHars;harN++)
927   {
928    fQ[selN][harN].Set(mQx[selN][harN],mQy[selN][harN]) ;  
929    for(subN=0;subN<AliFlowConstants::kSubs;subN++) { fQSub[subN][selN][harN].Set(mQxSub[subN][selN][harN],mQySub[subN][selN][harN]) ; }
930   }
931  }
932
933  fDone = kTRUE ;
934 }
935 //-----------------------------------------------------------------------