PWGCF/FLOW/Attic/AliFlowTrack.cxx

   1 //////////////////////////////////////////////////////////////////////
   2 //
   3 // $Id: AliFlowTrack.cxx 18618 2007-05-16 15:38:22Z snelling $
   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 //-------------------------------------------------------------