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