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