Adding comments (Laurent)
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTrack.cxx
CommitLineData
30a892e3 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
4e566f2f 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).
30a892e3 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"
92016a03 36
37#include "TMath.h"
30a892e3 38#include <iostream>
39using namespace std; //required for resolving the 'cout' symbol
40
92016a03 41ClassImp(AliFlowTrack)
42//-------------------------------------------------------------
d8b98e74 43AliFlowTrack::AliFlowTrack():
44fPhi(0.), fEta(0.), fPt(0.), fZFirstPoint(0.), fZLastPoint(0.), fChi2(0.), fTrackLength(0.), fMostLikelihoodPID(0), fPhiGlobal(0.), fEtaGlobal(0.), fPtGlobal(0.), fLabel(0) {
30a892e3 45 // Default constructor
46
4e566f2f 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++)
30a892e3 50 {
4e566f2f 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. ; }
30a892e3 54 }
55 ResetSelection() ;
56}
92016a03 57//-------------------------------------------------------------
d8b98e74 58AliFlowTrack::AliFlowTrack(const Char_t* name):
59fPhi(0.), fEta(0.), fPt(0.), fZFirstPoint(0.), fZLastPoint(0.), fChi2(0.), fTrackLength(0.), fMostLikelihoodPID(0), fPhiGlobal(0.), fEtaGlobal(0.), fPtGlobal(0.), fLabel(0) {
30a892e3 60 // TNamed constructor
61
62 SetName(name) ;
63 AliFlowTrack() ;
64}
92016a03 65//-------------------------------------------------------------
30a892e3 66AliFlowTrack::~AliFlowTrack()
67{
68 // default destructor (dummy)
69}
92016a03 70//-------------------------------------------------------------
30a892e3 71
92016a03 72//-------------------------------------------------------------
30a892e3 73Float_t AliFlowTrack::P() const
74{
92016a03 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.
30a892e3 77
92016a03 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
30a892e3 84 return momentum;
85}
92016a03 86//-------------------------------------------------------------
87Float_t AliFlowTrack::PGlobal() const
30a892e3 88{
92016a03 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) ; }
30a892e3 95
30a892e3 96 return momentum;
97}
92016a03 98//-------------------------------------------------------------
30a892e3 99Float_t AliFlowTrack::Mass() const
100{
101 // Returns the mass of the track, basing on its P.Id. hypotesis
102
92016a03 103 Float_t mass = 0.13957 ; // pion mass
30a892e3 104
92016a03 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; }
30a892e3 112
92016a03 113 return mass ;
30a892e3 114}
92016a03 115//-------------------------------------------------------------
30a892e3 116Float_t AliFlowTrack::InvMass() const
117{
118 // Returns the invariant mass of the track, calculated from T.O.F. and P_tot
119
92016a03 120 Float_t mass = -1 ; // if no calculation possible, it returns -1
30a892e3 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
92016a03 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() ; }
30a892e3 128 if(tof>0 && lenght>0)
92016a03 129 {
130 beta = lenght / (tof * c) ;
30a892e3 131 if(beta<1)
132 {
92016a03 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 ; }
30a892e3 139
92016a03 140 return mass ;
30a892e3 141}
92016a03 142//-------------------------------------------------------------
30a892e3 143Float_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 // -
92016a03 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)) ;
30a892e3 157
158 return rapidity ;
159}
92016a03 160//-------------------------------------------------------------
30a892e3 161Float_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 // -
92016a03 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)) ;
30a892e3 173
174 return rapidity ;
175}
92016a03 176//-------------------------------------------------------------
177const Char_t* AliFlowTrack::Pid() const
30a892e3 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
92016a03 182 const Char_t *name[] = {"e","mu","pi","k","pr","d"} ;
183 Int_t pdgCode = TMath::Abs(MostLikelihoodPID()) ;
184 Int_t charge = Charge() ;
30a892e3 185
92016a03 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" ; }
30a892e3 194
92016a03 195 if(charge>0) { pId += "+" ; }
196 else if(charge<0) { pId += "-" ; }
197 else { pId += "" ; }
30a892e3 198
92016a03 199 return pId.Data() ;
30a892e3 200}
92016a03 201//-------------------------------------------------------------
4e566f2f 202void AliFlowTrack::PidProbs(Float_t *pidN) const
30a892e3 203{
4e566f2f 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[]) .
30a892e3 207
208 Double_t sum = 0 ;
4e566f2f 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 ; }
30a892e3 213 }
214 else { cout << " ERROR - Empty Bayesian Vector !!! " << endl ; }
215}
92016a03 216//-------------------------------------------------------------
4e566f2f 217void 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//-------------------------------------------------------------
232Float_t AliFlowTrack::PidProb(Int_t nPid) const
30a892e3 233{
4e566f2f 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. ; }
30a892e3 238
4e566f2f 239 return (fCombRespFun[nPid] * AliFlowConstants::fgBayesian[nPid]) ;
30a892e3 240}
92016a03 241//-------------------------------------------------------------
4e566f2f 242Float_t AliFlowTrack::PidProbC(Int_t nPid) const
30a892e3 243{
4e566f2f 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 ...
30a892e3 246
4e566f2f 247 if(nPid > AliFlowConstants::kPid) { return 0. ; }
30a892e3 248
4e566f2f 249 Float_t customRespFun[AliFlowConstants::kPid] ;
250 GetCustomRespFun(customRespFun) ;
30a892e3 251
4e566f2f 252 return (customRespFun[nPid] * AliFlowConstants::fgBayesian[nPid]) ;
253}
92016a03 254//-------------------------------------------------------------
4e566f2f 255Float_t AliFlowTrack::MostLikelihoodRespFunc() const
30a892e3 256{
4e566f2f 257 // Returns the detector response function for the most probable P.id. hypothesis
258 // (Warning: THIS IS NOT THE BAYESIAN WEIGHT !)
30a892e3 259
92016a03 260 Int_t pdgCode = TMath::Abs(MostLikelihoodPID()) ;
4e566f2f 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] ; }
30a892e3 267 else { return 0. ; }
268}
92016a03 269//-------------------------------------------------------------
4e566f2f 270void 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//-------------------------------------------------------------
283void 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//-------------------------------------------------------------
296void 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//-------------------------------------------------------------
307void 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//-------------------------------------------------------------
30a892e3 325void 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}
92016a03 342//-------------------------------------------------------------
4e566f2f 343void 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//-------------------------------------------------------------
30a892e3 351Bool_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}
92016a03 362//-------------------------------------------------------------
363Float_t AliFlowTrack::Dca() const
30a892e3 364{
92016a03 365 // distance of closest approach (norm(dca[2]))
366
30a892e3 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}
92016a03 371//-------------------------------------------------------------
372Float_t AliFlowTrack::DcaError() const
30a892e3 373{
92016a03 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) ;
30a892e3 379}
92016a03 380//-------------------------------------------------------------
381void AliFlowTrack::DcaError3(Float_t err[3]) const
382{
383 // Dca Error (xy, z, cov(xy,z)) ...
30a892e3 384
92016a03 385 for(Int_t ii=0;ii<2;ii++) { err[ii] = TMath::Abs(fDcaError[ii]) ; }
386}
387//-------------------------------------------------------------
30a892e3 388
92016a03 389//-------------------------------------------------------------
390Bool_t AliFlowTrack::Select(Int_t harmonic,Int_t selection,Int_t subevent) const
391{
392 // Returns the selection flag for [harmonic] [selection] [sub-event]
30a892e3 393
92016a03 394 if((subevent == -1) || (subevent == fSubevent[harmonic][selection]))
395 {
396 return fSelection[harmonic][selection] ;
397 }
398 return kFALSE ;
399}
400//-------------------------------------------------------------
401void 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//-------------------------------------------------------------
417void 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)
30a892e3 421
92016a03 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//-------------------------------------------------------------
30a892e3 434void AliFlowTrack::SetConstrainable()
435{
4e566f2f 436 // fills the constrained parameters with the unconstrained ones,
437 // making the track a constrainable track.
438 // !!! TRICKY METHOD !!!
30a892e3 439
440 if(!IsConstrainable())
441 {
442 fPhi = fPhiGlobal ;
443 fEta = fEtaGlobal ;
444 fPt = fPtGlobal ;
445 }
446}
92016a03 447//-------------------------------------------------------------
30a892e3 448void AliFlowTrack::SetUnConstrainable()
449{
4e566f2f 450 // deletes the constrained parameters making the track an
451 // unconstrainable track.
452 // !!! TRICKY METHOD !!!
30a892e3 453
454 if(IsConstrainable())
455 {
456 fPhi = 0. ;
457 fEta = 0. ;
458 fPt = 0. ;
459 }
460}
92016a03 461//-------------------------------------------------------------