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