Improved version, kinetree included
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTrack.cxx
1 //////////////////////////////////////////////////////////////////////
2 //
3 // $Id$
4 //
5 // Author: Emanuele Simili
6 //
7 //////////////////////////////////////////////////////////////////////
8 //
9 //_____________________________________________________________
10 //
11 // Description: 
12 //         an array of AliFlowTrack is the core of the AliFlowEvent. 
13 // The object AliFlowTrack contains data members wich summarize the track 
14 // information most useful for flow study (such as Pt, eta, phi angle), some 
15 // detector signals (fitpoints, chi2 of the track fit, energy loss, time of 
16 // flight, ecc.), the p.id is stored as the detector response function, for 
17 // each detector (with the possibility to costomly combine them), and the 
18 // combined one (see bayesian P.Id. - chap.5 of the ALICE PPR). 
19 // This class is optimized for reaction plane calculation and sub-event 
20 // selection, through the appropriate methods in AliFlowEvent.  
21 // Two arrays of flags in the AliFlowTrack object (fSelection[][] and 
22 // fSubevent[][]) specify whether a track is included or not in a certain 
23 // Selection or Subevent. Those flags are switched on/off by the method 
24 // AliFlowEvent::SetSelection(), according to Eta, Pt, dca, constrainability  
25 // of the track, number of TPC hits, and p.id hypotesis (static member of the
26 // AliFlowEvent class (fEtaTpcCuts[][][], fPtTpcCuts[][][], fDcaGlobalCuts[], 
27 // fPid[], ecc.). 
28 //
29 // The AliFlowTrack Class is adapted from the original StFlowTrack,
30 // succesfully employed to study flow in the STAR experiment at RICH.
31 // Original Authors:                 Raimond Snellings & Art Poskanzer
32 //
33
34 #include "AliFlowTrack.h"
35 #include "AliFlowConstants.h"
36
37 #include "TMath.h"
38 #include <iostream>
39 using namespace std; //required for resolving the 'cout' symbol
40
41 ClassImp(AliFlowTrack)
42 //-------------------------------------------------------------
43 AliFlowTrack::AliFlowTrack()
44 {
45  // Default constructor
46  
47  fPhi = 0. ;
48  fEta = 0. ;
49  fPt  = 0. ;
50  fPhiGlobal = 0. ;
51  fEtaGlobal = 0. ;
52  fPtGlobal  = 0. ;
53  fChi2 = 0. ;
54  fZFirstPoint = 0. ;                  
55  fZLastPoint = 0. ;                           
56  fTrackLength = 0. ;
57  fMostLikelihoodPID = 0 ;                     
58  for(Int_t dd=0;dd<2;dd++)                      { fDcaSigned[dd] = 0. ; }                             
59  for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++) { fCombRespFun[ii]   = 0. ; }
60  for(Int_t det=0;det<4;det++)
61  {
62   fFitPts[det]  = 0  ; fMaxPts[det]   = 0  ; 
63   fFitChi2[det] = 0. ; fDedx[det]     = 0. ; fMom[det] = 0. ; 
64   for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++) { fRespFun[det][ii] = -1. ; }
65  }
66  ResetSelection() ;                           
67 }
68 //-------------------------------------------------------------
69 AliFlowTrack::AliFlowTrack(const Char_t* name) 
70 {
71  // TNamed constructor
72   
73  SetName(name) ;
74  AliFlowTrack() ;
75 }
76 //-------------------------------------------------------------
77 AliFlowTrack::~AliFlowTrack() 
78 {
79  // default destructor (dummy)
80 }
81 //-------------------------------------------------------------
82
83 //-------------------------------------------------------------
84 Float_t  AliFlowTrack::P() const 
85
86  // Returns the total momentum of the constrained track (calculated from Pt & Eta).
87  // If the track is not constrainable or eta is infinite, returns 0.
88  
89  if(!IsConstrainable()) { return 0. ; }               
90
91  Float_t momentum = 0. ; 
92  Float_t conv = 1 - ( TMath::TanH(Eta()) * TMath::TanH(Eta()) ) ;
93  if(conv>0) { momentum = Pt() / TMath::Sqrt(conv) ; } 
94
95  return momentum; 
96 }
97 //-------------------------------------------------------------
98 Float_t  AliFlowTrack::PGlobal() const 
99
100  // Returns the total momentum of the unconstrained track (calculated from PtGlobal & EtaGlobal). 
101  // If etaGlobal is infinite, returns 0.
102
103  Float_t momentum = 0. ; 
104  Float_t conv = 1 - ( TMath::TanH(EtaGlobal()) * TMath::TanH(EtaGlobal()) ) ;
105  if(conv>0) { momentum = PtGlobal() / TMath::Sqrt(conv) ; }
106
107  return momentum; 
108 }
109 //-------------------------------------------------------------
110 Float_t  AliFlowTrack::Mass() const 
111
112  // Returns the mass of the track, basing on its P.Id. hypotesis
113
114  Float_t mass = 0.13957 ;  // pion mass 
115
116  if(MostLikelihoodPID() == 0)                         { mass = -1.    ; }
117  else if(TMath::Abs(MostLikelihoodPID()) == 11)       { mass = 0.00051; }
118  else if(TMath::Abs(MostLikelihoodPID()) == 13)       { mass = 0.10566; }
119  else if(TMath::Abs(MostLikelihoodPID()) == 211)      { mass = 0.49368; }
120  else if(TMath::Abs(MostLikelihoodPID()) == 321)      { mass = 0.13957; }
121  else if(TMath::Abs(MostLikelihoodPID()) == 2212)     { mass = 0.93827; }
122  else if(TMath::Abs(MostLikelihoodPID()) == 10010020) { mass = 1.87505; }
123
124  return mass ;
125 }
126 //-------------------------------------------------------------
127 Float_t  AliFlowTrack::InvMass() const 
128
129  // Returns the invariant mass of the track, calculated from T.O.F. and P_tot
130
131  Float_t mass = -1 ;                // if no calculation possible, it returns -1
132  Float_t c = TMath::Ccgs()/1e+12 ;  // = 0.02998 ;  --> speed of light in cm/psec
133  Float_t tof = TofTOF() ;           // in pico-seconds
134  Float_t lenght = TrackLength() ;   // in cm
135  Float_t ptot = 0 ;                 // constrained parameters are used if there, otherwise unconstrained
136  Float_t beta, gamma ;
137  if(IsConstrainable()) { ptot = P() ; }
138  else                  { ptot = PGlobal() ; }
139  if(tof>0 && lenght>0)
140  {
141   beta = lenght / (tof * c) ;
142   if(beta<1) 
143   {
144    gamma = 1 / TMath::Sqrt( 1 - (beta * beta) ) ;
145    mass = TMath::Abs( ptot / (beta * gamma) ) ;  
146   }
147   //cout << " (yes) Mass = " << mass << " , t/l = " << (tof/lenght) << " , beta = " << beta << " , gamma = " << gamma << endl ; 
148  }   
149  //else { cout << " (no)  TOF = " << tof << "   , Lenght = " << lenght << endl ; }
150
151  return mass ;
152 }
153 //-------------------------------------------------------------
154 Float_t  AliFlowTrack::Y() const 
155 {
156  // Rapidity of the constrained track.
157  // If track is not constrainable, the unconstrained rapidity is returned.
158  
159  if(!IsConstrainable()) { return YGlobal() ; }
160  if(TMath::Abs((Float_t)P()) == TMath::Abs((Float_t)Pt()))      { return 0. ; }
161  else if(TMath::Abs((Float_t)P()) < TMath::Abs((Float_t)Pt()))  { cout << "track: " << GetName() << "has  Pt() > P() !!!" << endl ; }
162  // -
163  Float_t mass = Mass() ; 
164  Double_t pz = TMath::Sqrt(P()*P() - Pt()*Pt()); 
165  if(Eta()<0) { pz = -pz ; }
166  Double_t e = TMath::Sqrt(P()*P() + mass*mass) ;
167  Float_t rapidity = 0.5*TMath::Log((e + pz)/(e - pz)) ;
168  
169  return rapidity ;
170 }
171 //-------------------------------------------------------------
172 Float_t  AliFlowTrack::YGlobal() const 
173 {
174  // Rapidity of the unconstrained track
175
176  if(TMath::Abs((Float_t)PGlobal()) == TMath::Abs((Float_t)PtGlobal()))     { return 0. ; }
177  else if(TMath::Abs((Float_t)PGlobal()) < TMath::Abs((Float_t)PtGlobal())) { cout << "track: " << GetName() << "has  Pt() > P() !!!" << endl ; }
178  // -
179  Float_t mass = Mass() ; 
180  Double_t pz = TMath::Sqrt(PGlobal()*PGlobal() - PtGlobal()*PtGlobal()); 
181  if(EtaGlobal()<0) { pz = -pz ; }
182  Double_t e = TMath::Sqrt(PGlobal()*PGlobal() + mass*mass) ;
183  Float_t rapidity = 0.5*TMath::Log((e + pz)/(e - pz)) ;
184  
185  return rapidity ;
186 }
187 //-------------------------------------------------------------
188 const Char_t* AliFlowTrack::Pid() const
189 {
190  // Returns the P.Id. in characters (e,mu,pi,k,p,d) basing on the stored pdg code 
191  // and its sign (both stored in MostLikelihoodPID() ) .
192  
193  const Char_t *name[] = {"e","mu","pi","k","pr","d"} ;
194  Int_t pdgCode = TMath::Abs(MostLikelihoodPID()) ;
195  Int_t charge = Charge() ;
196
197  TString pId = "" ;
198  if(pdgCode == 11)            { pId = name[0] ; }
199  else if(pdgCode == 13)       { pId = name[1] ; }
200  else if(pdgCode == 211)      { pId = name[2] ; }
201  else if(pdgCode == 321)      { pId = name[3] ; }
202  else if(pdgCode == 2212)     { pId = name[4] ; }
203  else if(pdgCode == 10010020) { pId = name[5] ; }
204  else                         { pId = "0" ; }
205
206  if(charge>0)                 { pId += "+" ; } 
207  else if(charge<0)            { pId += "-" ; }
208  else                         { pId += "" ; }
209  
210  return pId.Data() ; 
211 }
212 //-------------------------------------------------------------
213 void AliFlowTrack::PidProbs(Float_t *pidN) const 
214
215  // Returns the normalized probability (the "bayesian weights") for the given track to be [e,mu,pi,k,p,d] .  
216  // The COMBINED detector response function is scaled by the a priori probabilities 
217  // (the normalised particle abundances is stored in AliFlowConstants::fgBayesian[]) .
218
219  Double_t sum = 0 ; 
220  for(Int_t n=0;n<AliFlowConstants::kPid;n++) { pidN[n] = PidProb(n) ; sum += pidN[n] ; }
221  if(sum) 
222  {  
223   for(Int_t n=0;n<AliFlowConstants::kPid;n++)  { pidN[n] /=  sum ; }
224  }
225  else { cout << " ERROR - Empty Bayesian Vector !!! " << endl ; }
226
227 //-------------------------------------------------------------
228 void AliFlowTrack::PidProbsC(Float_t *pidN) const 
229
230  // Returns the normalized probability (the "bayesian weights") for the given track to be [e,mu,pi,k,p,d] .  
231  // The CUSTOM detector response function (see AliFlowTrack) is scaled by the a priori probabilities 
232  // (the normalised particle abundances is stored in AliFlowConstants::fgBayesian[]) .
233  
234  Double_t sum = 0 ; 
235  for(Int_t n=0;n<AliFlowConstants::kPid;n++) { pidN[n] = PidProbC(n) ; sum += pidN[n] ; }
236  if(sum) 
237  {  
238   for(Int_t n=0;n<AliFlowConstants::kPid;n++)  { pidN[n] /=  sum ; }
239  }
240  else { cout << " ERROR - Empty Bayesian Vector !!! " << endl ; }
241
242 //-------------------------------------------------------------
243 Float_t  AliFlowTrack::PidProb(Int_t nPid) const
244 {
245  // Returns the bayesian weight of the track to be [nPid = e,mu,pi,k,pi,d].
246  // The detector response function in use is the combined one (from the ESD)
247
248  if(nPid > AliFlowConstants::kPid) { return 0. ; } 
249
250  return (fCombRespFun[nPid] * AliFlowConstants::fgBayesian[nPid]) ;
251 }
252 //-------------------------------------------------------------
253 Float_t  AliFlowTrack::PidProbC(Int_t nPid) const
254 {
255  // Returns the bayesian weight of the track to be [nPid = e,mu,pi,k,pi,d].
256  // The detector response function in use is the custom one ...
257
258  if(nPid > AliFlowConstants::kPid) { return 0. ; } 
259
260  Float_t  customRespFun[AliFlowConstants::kPid] ;
261  GetCustomRespFun(customRespFun) ;
262
263  return (customRespFun[nPid] * AliFlowConstants::fgBayesian[nPid]) ;
264 }
265 //-------------------------------------------------------------
266 Float_t AliFlowTrack::MostLikelihoodRespFunc() const 
267
268  // Returns the detector response function for the most probable P.id. hypothesis
269  // (Warning: THIS IS NOT THE BAYESIAN WEIGHT !) 
270
271  Int_t pdgCode = TMath::Abs(MostLikelihoodPID()) ;
272  if(pdgCode == 11)            { return fCombRespFun[0] ; }
273  else if(pdgCode == 13)       { return fCombRespFun[1] ; }
274  else if(pdgCode == 211)      { return fCombRespFun[2] ; }
275  else if(pdgCode == 321)      { return fCombRespFun[3] ; }
276  else if(pdgCode == 2212)     { return fCombRespFun[4] ; }
277  else if(pdgCode == 10010020) { return fCombRespFun[5] ; }
278  else { return 0. ; }
279
280 //-------------------------------------------------------------
281 void AliFlowTrack::SetRespFun(Int_t det, Float_t *r)
282 {
283  // This function fills "AliFlowConstants::kPid" PID weights 
284  // (detector response functions) of the track .
285  // The method is private, cause it is called by the public 
286  // methods: SetRespFunITS, SetRespFunTPC, ...
287
288  for(Int_t i=0;i<AliFlowConstants::kPid;i++)
289  {
290   fRespFun[det][i] = r[i] ;
291  }
292 }
293 //-------------------------------------------------------------
294 void AliFlowTrack::GetRespFun(Int_t det, Float_t *r) const
295 {
296  // This function returns the response functions of the 
297  // detector [det] for the P.Id. of the track .
298  // The method is private, and it is called by the public 
299  // methods: GetRespFunITS, GetRespFunTPC, ...
300
301  for(Int_t i=0;i<AliFlowConstants::kPid;i++)
302  {
303   r[i] = fRespFun[det][i] ;
304  }
305 }
306 //-------------------------------------------------------------
307 void AliFlowTrack::GetCombinedRespFun(Float_t *r) const 
308 {
309  // This function returns the combined response functions for  
310  // [e,mu,pi,k,pi,d] as it is stored in the ESD .
311  
312  for(Int_t i=0;i<AliFlowConstants::kPid;i++)
313  {
314   r[i] = fCombRespFun[i] ;
315  }
316 }
317 //------------------------------------------------------------- 
318 void AliFlowTrack::GetCustomRespFun(Float_t *r) const
319 {
320  // This function returns a combined response functions as setted 
321  // by the user ... at the moment just the sum of single responses ...  
322  // for the track to be [e,mu,pi,k,pi,d] .
323  
324  Int_t sum ;
325  for(Int_t i=0;i<AliFlowConstants::kPid;i++)
326  {
327   r[i] = 0. ; sum = 0 ;
328   for(Int_t det=0;det<4;det++) 
329   { 
330    if(fRespFun[det][i]>0) { r[i] += fRespFun[det][i] ; sum += 1 ; } 
331   }
332   if(sum) { r[i] /= sum ; }
333  }
334 }
335 //-------------------------------------------------------------
336 void AliFlowTrack::SetPid(const Char_t* pid)            
337
338  // Sets the P.Id. hypotesis of the track from a String imput (that should be 
339  // ("e","mu","pi","k","pr","d"). Sign is allowed as well. The string itself 
340  // is not stored, what is stored is the signed PDG code.
341  
342  if(strstr(pid,"e"))       { fMostLikelihoodPID = 11    ; }
343  else if(strstr(pid,"mu")) { fMostLikelihoodPID = 13    ; }
344  else if(strstr(pid,"pi")) { fMostLikelihoodPID = 211   ; }
345  else if(strstr(pid,"k"))  { fMostLikelihoodPID = 321   ; }
346  else if(strstr(pid,"pr")) { fMostLikelihoodPID = 2212  ; }
347  else if(strstr(pid,"d"))  { fMostLikelihoodPID = 10010020 ; }
348  else { fMostLikelihoodPID = 0 ; cout << "AliFlowTrack - !BAD IMPUT!" << endl ; }
349  
350  if(strchr(pid,'+'))       { fMostLikelihoodPID = TMath::Abs(fMostLikelihoodPID) * 1 ; }
351  else if(strchr(pid,'-'))  { fMostLikelihoodPID = TMath::Abs(fMostLikelihoodPID) * -1 ; } 
352 }
353 //-------------------------------------------------------------
354 void AliFlowTrack::SetPid(Int_t pdgCode)                
355
356  // Sets the P.Id. hypotesis of the track from the PDG code. 
357  // Sign can be given as well. 
358  
359  fMostLikelihoodPID = pdgCode ;
360 }
361 //-------------------------------------------------------------
362 Bool_t AliFlowTrack::IsConstrainable()              const 
363
364  // Returns kTRUE if the track is constrainable. 
365  // If the track is constrainable -> the constrained parameters are stored (and the
366  // unconstrained ones as well in ...Global). 
367  // Otherwise (if track is not constrainable) -> just the unconstrained parameters 
368  // are stored, and returned instead of the constrained ones .
369  
370  if(fPt==0 && fEta==0) { return kFALSE ; }
371  else                  { return kTRUE ; }
372
373 //-------------------------------------------------------------
374 Float_t AliFlowTrack::Dca() const 
375 {
376  // distance of closest approach (norm(dca[2]))
377
378  Double_t mag = 0 ;
379  for(Int_t ii=0;ii<2;ii++) { mag += (fDcaSigned[ii]*fDcaSigned[ii]) ; } 
380  return TMath::Sqrt(mag) ; 
381 }     
382 //-------------------------------------------------------------
383 Float_t AliFlowTrack::DcaError() const 
384 {
385  // error on the distance of closest approach
386
387  Double_t err = 0 ;
388  for(Int_t ii=0;ii<2;ii++) { err = (fDcaError[ii]*fDcaError[ii]) ; } 
389  return TMath::Sqrt(err) ; 
390 }     
391 //-------------------------------------------------------------
392 void AliFlowTrack::DcaError3(Float_t err[3]) const 
393 {
394  // Dca Error (xy, z, cov(xy,z)) ...
395
396  for(Int_t ii=0;ii<2;ii++) { err[ii] = TMath::Abs(fDcaError[ii]) ; }            
397 }     
398 //-------------------------------------------------------------
399
400 //-------------------------------------------------------------
401 Bool_t AliFlowTrack::Select(Int_t harmonic,Int_t selection,Int_t subevent) const 
402 {
403  // Returns the selection flag for [harmonic] [selection] [sub-event]
404
405  if((subevent == -1) || (subevent == fSubevent[harmonic][selection])) 
406  {
407   return fSelection[harmonic][selection] ; 
408  } 
409  return kFALSE ;        
410 }
411 //-------------------------------------------------------------
412 void AliFlowTrack::PrintSelection() const
413 {
414  // Prints a short string of 0 & 1 to visualize the selections' flags
415  // [har1][sel0],[har1][sel1],[har1][sel2],...,[har2][sel0],...
416  
417  for(Int_t i=0;i<AliFlowConstants::kHars;i++)
418  {
419   for(Int_t j=0;j<AliFlowConstants::kSels;j++)
420   {
421    if(Select(i,j)) { cout << 1 ; }
422    else            { cout << 0 ; }
423   }
424  }
425  cout << endl ;
426 }
427 //-------------------------------------------------------------
428 void AliFlowTrack::ResetSelection() 
429 {
430  // Re-sets all the selection/sub-event flags to 0 
431  // (track won't be used in any R.P. calculation)
432
433  for(Int_t ii=0;ii<AliFlowConstants::kHars;ii++)
434  {
435   for(Int_t jj=0;jj<AliFlowConstants::kSels;jj++)
436   {
437    fSelection[ii][jj] = kFALSE ; 
438    fSubevent[ii][jj] = -1 ; 
439   }
440  }
441 }
442 //-------------------------------------------------------------
443
444 //-------------------------------------------------------------
445 void AliFlowTrack::SetConstrainable()                         
446
447  // fills the constrained parameters with the unconstrained ones,
448  // making the track a constrainable track.
449  //                                       !!! TRICKY METHOD !!!
450
451  if(!IsConstrainable()) 
452  { 
453   fPhi = fPhiGlobal ; 
454   fEta = fEtaGlobal ; 
455   fPt  = fPtGlobal ; 
456  } 
457 }
458 //-------------------------------------------------------------
459 void AliFlowTrack::SetUnConstrainable()                       
460
461  // deletes the constrained parameters making the track an 
462  // unconstrainable track. 
463  //                                       !!! TRICKY METHOD !!!
464
465  if(IsConstrainable()) 
466  { 
467   fPhi = 0. ; 
468   fEta = 0. ; 
469   fPt  = 0. ; 
470  } 
471 }
472 //-------------------------------------------------------------