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