STEER/ESD/AliKFParticle.cxx

   1 //----------------------------------------------------------------------------
   2 // Implementation of the AliKFParticle class
   3 // .
   4 // @author  S.Gorbunov, I.Kisel
   5 // @version 1.0
   6 // @since   13.05.07
   7 //
   8 // Class to reconstruct and store the decayed particle parameters.
   9 // The method is described in CBM-SOFT note 2007-003,
  10 // ``Reconstruction of decayed particles based on the Kalman filter'',
  11 // http://www.gsi.de/documents/DOC-2007-May-14-1.pdf
  12 //
  13 // This class is ALICE interface to general mathematics in AliKFParticleCore
  14 //
  15 //  -= Copyright &copy ALICE HLT Group =-
  16 //____________________________________________________________________________
  17
  18
  19 #include "AliKFParticle.h"
  20 #include "TDatabasePDG.h"
  21 #include "TParticlePDG.h"
  22 #include "AliVTrack.h"
  23 #include "AliVVertex.h"
  24
  25 #include "AliExternalTrackParam.h"
  26
  27 ClassImp(AliKFParticle)
  28
  29 Double_t AliKFParticle::fgBz = -5.;  //* Bz compoment of the magnetic field
  30
  31 AliKFParticle::AliKFParticle( const AliKFParticle &d1, const AliKFParticle &d2, Bool_t gamma )
  32 {
  33   if (!gamma) {
  34     AliKFParticle mother;
  35     mother+= d1;
  36     mother+= d2;
  37     *this = mother;
  38   } else
  39     ConstructGamma(d1, d2);
  40 }
  41
  42 void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
  43 {
  44   // Constructor from "cartesian" track, PID hypothesis should be provided
  45   //
  46   // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
  47   // Cov [21] = lower-triangular part of the covariance matrix:
  48   //
  49   //                (  0  .  .  .  .  . )
  50   //                (  1  2  .  .  .  . )
  51   //  Cov. matrix = (  3  4  5  .  .  . ) - numbering of covariance elements in Cov[]
  52   //                (  6  7  8  9  .  . )
  53   //                ( 10 11 12 13 14  . )
  54   //                ( 15 16 17 18 19 20 )
  55   Double_t C[21];
  56   for( int i=0; i<21; i++ ) C[i] = Cov[i];
  57
  58   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
  59   Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
  60
  61   AliKFParticleBase::Initialize( Param, C, Charge, mass );
  62 }
  63
  64 void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Double_t Mass )
  65 {
  66   // Constructor from "cartesian" track, PID hypothesis should be provided
  67   //
  68   // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
  69   // Cov [21] = lower-triangular part of the covariance matrix:
  70   //
  71   //                (  0  .  .  .  .  . )
  72   //                (  1  2  .  .  .  . )
  73   //  Cov. matrix = (  3  4  5  .  .  . ) - numbering of covariance elements in Cov[]
  74   //                (  6  7  8  9  .  . )
  75   //                ( 10 11 12 13 14  . )
  76   //                ( 15 16 17 18 19 20 )
  77   Double_t C[21];
  78   for( int i=0; i<21; i++ ) C[i] = Cov[i];
  79
  80   AliKFParticleBase::Initialize( Param, C, Charge, Mass );
  81 }
  82
  83 AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
  84 {
  85   // Constructor from ALICE track, PID hypothesis should be provided
  86
  87   track.GetXYZ(fP);
  88   track.PxPyPz(fP+3);
  89   fQ = track.Charge();
  90   track.GetCovarianceXYZPxPyPz( fC );
  91   Create(fP,fC,fQ,PID);
  92 }
  93
  94 AliKFParticle::AliKFParticle( const AliExternalTrackParam &track, Double_t Mass, Int_t Charge )
  95 {
  96   // Constructor from ALICE track, Mass and Charge hypothesis should be provided
  97
  98   track.GetXYZ(fP);
  99   track.PxPyPz(fP+3);
 100   fQ = track.Charge()*TMath::Abs(Charge);
 101   fP[3] *= TMath::Abs(Charge);
 102   fP[4] *= TMath::Abs(Charge);
 103   fP[5] *= TMath::Abs(Charge);
 104
 105   Double_t pt=1./TMath::Abs(track.GetParameter()[4]) * TMath::Abs(Charge);
 106   Double_t cs=TMath::Cos(track.GetAlpha()), sn=TMath::Sin(track.GetAlpha());
 107   Double_t r=TMath::Sqrt((1.-track.GetParameter()[2])*(1.+track.GetParameter()[2]));
 108
 109   Double_t m00=-sn, m10=cs;
 110   Double_t m23=-pt*(sn + track.GetParameter()[2]*cs/r), m43=-pt*pt*(r*cs - track.GetParameter()[2]*sn);
 111   Double_t m24= pt*(cs - track.GetParameter()[2]*sn/r), m44=-pt*pt*(r*sn + track.GetParameter()[2]*cs);
 112   Double_t m35=pt, m45=-pt*pt*track.GetParameter()[3];
 113
 114   m43*=track.GetSign();
 115   m44*=track.GetSign();
 116   m45*=track.GetSign();
 117
 118   const Double_t *cTr = track.GetCovariance();
 119
 120   fC[0 ] = cTr[0]*m00*m00;
 121   fC[1 ] = cTr[0]*m00*m10;
 122   fC[2 ] = cTr[0]*m10*m10;
 123   fC[3 ] = cTr[1]*m00;
 124   fC[4 ] = cTr[1]*m10;
 125   fC[5 ] = cTr[2];
 126   fC[6 ] = m00*(cTr[3]*m23 + cTr[10]*m43);
 127   fC[7 ] = m10*(cTr[3]*m23 + cTr[10]*m43);
 128   fC[8 ] = cTr[4]*m23 + cTr[11]*m43;
 129   fC[9 ] = m23*(cTr[5]*m23 + cTr[12]*m43)  +  m43*(cTr[12]*m23 + cTr[14]*m43);
 130   fC[10] = m00*(cTr[3]*m24 + cTr[10]*m44);
 131   fC[11] = m10*(cTr[3]*m24 + cTr[10]*m44);
 132   fC[12] = cTr[4]*m24 + cTr[11]*m44;
 133   fC[13] = m23*(cTr[5]*m24 + cTr[12]*m44)  +  m43*(cTr[12]*m24 + cTr[14]*m44);
 134   fC[14] = m24*(cTr[5]*m24 + cTr[12]*m44)  +  m44*(cTr[12]*m24 + cTr[14]*m44);
 135   fC[15] = m00*(cTr[6]*m35 + cTr[10]*m45);
 136   fC[16] = m10*(cTr[6]*m35 + cTr[10]*m45);
 137   fC[17] = cTr[7]*m35 + cTr[11]*m45;
 138   fC[18] = m23*(cTr[8]*m35 + cTr[12]*m45)  +  m43*(cTr[13]*m35 + cTr[14]*m45);
 139   fC[19] = m24*(cTr[8]*m35 + cTr[12]*m45)  +  m44*(cTr[13]*m35 + cTr[14]*m45);
 140   fC[20] = m35*(cTr[9]*m35 + cTr[13]*m45)  +  m45*(cTr[13]*m35 + cTr[14]*m45);
 141
 142   Create(fP,fC,fQ,Mass);
 143 }
 144
 145 AliKFParticle::AliKFParticle( const AliVVertex &vertex )
 146 {
 147   // Constructor from ALICE vertex
 148
 149   vertex.GetXYZ( fP );
 150   vertex.GetCovarianceMatrix( fC );
 151   fChi2 = vertex.GetChi2();
 152   fNDF = 2*vertex.GetNContributors() - 3;
 153   fQ = 0;
 154   fAtProductionVertex = 0;
 155   fIsLinearized = 0;
 156   fSFromDecay = 0;
 157 }
 158
 159 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
 160 {
 161   // Conversion to AliExternalTrackParam parameterization
 162
 163   Double_t cosA = p.GetPx(), sinA = p.GetPy();
 164   Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
 165   Double_t pti = 0;
 166   if( pt<1.e-4 ){
 167     cosA = 1;
 168     sinA = 0;
 169   } else {
 170     pti = 1./pt;
 171     cosA*=pti;
 172     sinA*=pti;
 173   }
 174   Alpha = TMath::ATan2(sinA,cosA);
 175   X   = p.GetX()*cosA + p.GetY()*sinA;
 176   P[0]= p.GetY()*cosA - p.GetX()*sinA;
 177   P[1]= p.GetZ();
 178   P[2]= 0;
 179   P[3]= p.GetPz()*pti;
 180   P[4]= p.GetQ()*pti;
 181 }
 182
 183 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
 184 {
 185   //* Calculate DCA distance from vertex (transverse impact parameter) in XY
 186   //* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix
 187
 188   Bool_t ret = 0;
 189
 190   Double_t mP[8];
 191   Double_t mC[36];
 192
 193   Transport( GetDStoPoint(vtx), mP, mC );
 194
 195   Double_t dx = mP[0] - vtx[0];
 196   Double_t dy = mP[1] - vtx[1];
 197   Double_t px = mP[3];
 198   Double_t py = mP[4];
 199   Double_t pt = TMath::Sqrt(px*px + py*py);
 200   Double_t ex=0, ey=0;
 201   if( pt<1.e-4 ){
 202     ret = 1;
 203     pt = 1.;
 204     val = 1.e4;
 205   } else{
 206     ex = px/pt;
 207     ey = py/pt;
 208     val = dy*ex - dx*ey;
 209   }
 210
 211   Double_t h0 = -ey;
 212   Double_t h1 = ex;
 213   Double_t h3 = (dy*ey + dx*ex)*ey/pt;
 214   Double_t h4 = -(dy*ey + dx*ex)*ex/pt;
 215
 216   err =
 217     h0*(h0*GetCovariance(0,0) + h1*GetCovariance(0,1) + h3*GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
 218     h1*(h0*GetCovariance(1,0) + h1*GetCovariance(1,1) + h3*GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
 219     h3*(h0*GetCovariance(3,0) + h1*GetCovariance(3,1) + h3*GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
 220     h4*(h0*GetCovariance(4,0) + h1*GetCovariance(4,1) + h3*GetCovariance(4,3) + h4*GetCovariance(4,4) );
 221
 222   if( Cv ){
 223     err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] );
 224   }
 225
 226   err = TMath::Sqrt(TMath::Abs(err));
 227
 228   return ret;
 229 }
 230
 231 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
 232 {
 233   return GetDistanceFromVertexXY( vtx, 0, val, err );
 234 }
 235
 236
 237 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const
 238 {
 239   //* Calculate distance from vertex [cm] in XY-plane
 240
 241   return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
 242 }
 243
 244 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx, Double_t &val, Double_t &err ) const
 245 {
 246   //* Calculate distance from vertex [cm] in XY-plane
 247
 248   return GetDistanceFromVertexXY( AliKFParticle(Vtx), val, err );
 249 }
 250
 251 Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
 252 {
 253   //* Calculate distance from vertex [cm] in XY-plane
 254   Double_t val, err;
 255   GetDistanceFromVertexXY( vtx, 0, val, err );
 256   return val;
 257 }
 258
 259 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const
 260 {
 261   //* Calculate distance from vertex [cm] in XY-plane
 262
 263   return GetDistanceFromVertexXY( Vtx.fP );
 264 }
 265
 266 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const
 267 {
 268   //* Calculate distance from vertex [cm] in XY-plane
 269
 270   return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
 271 }
 272
 273 Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const
 274 {
 275   //* Calculate distance to other particle [cm]
 276
 277   Double_t dS, dS1;
 278   GetDStoParticleXY( p, dS, dS1 );
 279   Double_t mP[8], mC[36], mP1[8], mC1[36];
 280   Transport( dS, mP, mC );
 281   p.Transport( dS1, mP1, mC1 );
 282   Double_t dx = mP[0]-mP1[0];
 283   Double_t dy = mP[1]-mP1[1];
 284   return TMath::Sqrt(dx*dx+dy*dy);
 285 }
 286
 287 Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const
 288 {
 289   //* Calculate sqrt(Chi2/ndf) deviation from other particle
 290
 291   Double_t dS, dS1;
 292   GetDStoParticleXY( p, dS, dS1 );
 293   Double_t mP1[8], mC1[36];
 294   p.Transport( dS1, mP1, mC1 );
 295
 296   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
 297
 298   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
 299                                         (mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );
 300
 301   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };
 302
 303   mC1[0] +=h[0]*h[0];
 304   mC1[1] +=h[1]*h[0];
 305   mC1[2] +=h[1]*h[1];
 306
 307   return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
 308 }
 309
 310
 311 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const
 312 {
 313   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 314   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 315
 316   Double_t val, err;
 317   Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
 318   if( problem || err<1.e-20 ) return 1.e4;
 319   else return val/err;
 320 }
 321
 322
 323 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
 324 {
 325   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 326   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 327
 328   return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
 329 }
 330
 331 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const
 332 {
 333   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 334   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 335
 336   AliKFParticle v(Vtx);
 337   return GetDeviationFromVertexXY( v.fP, v.fC );
 338 }
 339
 340 Double_t AliKFParticle::GetAngle  ( const AliKFParticle &p ) const
 341 {
 342   //* Calculate the opening angle between two particles
 343
 344   Double_t dS, dS1;
 345   GetDStoParticle( p, dS, dS1 );
 346   Double_t mP[8], mC[36], mP1[8], mC1[36];
 347   Transport( dS, mP, mC );
 348   p.Transport( dS1, mP1, mC1 );
 349   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
 350   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
 351   n*=n1;
 352   Double_t a = 0;
 353   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
 354   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 355   else a = (a>=0) ?0 :TMath::Pi();
 356   return a;
 357 }
 358
 359 Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const
 360 {
 361   //* Calculate the opening angle between two particles in XY plane
 362
 363   Double_t dS, dS1;
 364   GetDStoParticleXY( p, dS, dS1 );
 365   Double_t mP[8], mC[36], mP1[8], mC1[36];
 366   Transport( dS, mP, mC );
 367   p.Transport( dS1, mP1, mC1 );
 368   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 369   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
 370   n*=n1;
 371   Double_t a = 0;
 372   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
 373   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 374   else a = (a>=0) ?0 :TMath::Pi();
 375   return a;
 376 }
 377
 378 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const
 379 {
 380   //* Calculate the opening angle between two particles in RZ plane
 381
 382   Double_t dS, dS1;
 383   GetDStoParticle( p, dS, dS1 );
 384   Double_t mP[8], mC[36], mP1[8], mC1[36];
 385   Transport( dS, mP, mC );
 386   p.Transport( dS1, mP1, mC1 );
 387   Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 388   Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
 389   Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
 390   Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
 391   n*=n1;
 392   Double_t a = 0;
 393   if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n;
 394   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 395   else a = (a>=0) ?0 :TMath::Pi();
 396   return a;
 397 }
 398
 399
 400 /*
 401
 402 #include "AliExternalTrackParam.h"
 403
 404 void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p,
 405                                            Double_t &DS, Double_t &DS1 )
 406   const
 407 {
 408   DS = DS1 = 0;
 409   Double_t x1, a1, x2, a2;
 410   Double_t par1[5], par2[5], cov[15];
 411   for(int i=0; i<15; i++) cov[i] = 0;
 412   cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;
 413
 414   GetExternalTrackParam( *this, x1, a1, par1 );
 415   GetExternalTrackParam( p, x2, a2, par2 );
 416
 417   AliExternalTrackParam t1(x1,a1, par1, cov);
 418   AliExternalTrackParam t2(x2,a2, par2, cov);
 419
 420   Double_t xe1=0, xe2=0;
 421   t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
 422   t1.PropagateTo( xe1, -GetFieldAlice() );
 423   t2.PropagateTo( xe2, -GetFieldAlice() );
 424
 425   Double_t xyz1[3], xyz2[3];
 426   t1.GetXYZ( xyz1 );
 427   t2.GetXYZ( xyz2 );
 428
 429   DS = GetDStoPoint( xyz1 );
 430   DS1 = p.GetDStoPoint( xyz2 );
 431
 432   return;
 433 }
 434 */
 435
 436   // * Pseudo Proper Time of decay = (r*pt) / |pt| * M/|pt|
 437 Double_t AliKFParticle::GetPseudoProperDecayTime( const AliKFParticle &pV, const Double_t& mass, Double_t* timeErr2 ) const
 438 { // TODO optimize with respect to time and stability
 439   const Double_t ipt2 = 1/( Px()*Px() + Py()*Py() );
 440   const Double_t mipt2 = mass*ipt2;
 441   const Double_t dx = X() - pV.X();
 442   const Double_t dy = Y() - pV.Y();
 443
 444   if ( timeErr2 ) {
 445       // -- calculate error = sigma(f(r)) = f'Cf'
 446       // r = {x,y,px,py,x_pV,y_pV}
 447       // df/dr = { px*m/pt^2,
 448       //     py*m/pt^2,
 449       //    ( x - x_pV )*m*(1/pt^2 - 2(px/pt^2)^2),
 450       //    ( y - y_pV )*m*(1/pt^2 - 2(py/pt^2)^2),
 451       //     -px*m/pt^2,
 452       //     -py*m/pt^2 }
 453     const Double_t f0 = Px()*mipt2;
 454     const Double_t f1 = Py()*mipt2;
 455     const Double_t mipt2derivative = mipt2*(1-2*Px()*Px()*ipt2);
 456     const Double_t f2 = dx*mipt2derivative;
 457     const Double_t f3 = -dy*mipt2derivative;
 458     const Double_t f4 = -f0;
 459     const Double_t f5 = -f1;
 460
 461     const Double_t& mC00 =    GetCovariance(0,0);
 462     const Double_t& mC10 =    GetCovariance(0,1);
 463     const Double_t& mC11 =    GetCovariance(1,1);
 464     const Double_t& mC20 =    GetCovariance(3,0);
 465     const Double_t& mC21 =    GetCovariance(3,1);
 466     const Double_t& mC22 =    GetCovariance(3,3);
 467     const Double_t& mC30 =    GetCovariance(4,0);
 468     const Double_t& mC31 =    GetCovariance(4,1);
 469     const Double_t& mC32 =    GetCovariance(4,3);
 470     const Double_t& mC33 =    GetCovariance(4,4);
 471     const Double_t& mC44 = pV.GetCovariance(0,0);
 472     const Double_t& mC54 = pV.GetCovariance(1,0);
 473     const Double_t& mC55 = pV.GetCovariance(1,1);
 474
 475     *timeErr2 =
 476       f5*mC55*f5 +
 477       f5*mC54*f4 +
 478       f4*mC44*f4 +
 479       f3*mC33*f3 +
 480       f3*mC32*f2 +
 481       f3*mC31*f1 +
 482       f3*mC30*f0 +
 483       f2*mC22*f2 +
 484       f2*mC21*f1 +
 485       f2*mC20*f0 +
 486       f1*mC11*f1 +
 487       f1*mC10*f0 +
 488       f0*mC00*f0;
 489   }
 490   return ( dx*Px() + dy*Py() )*mipt2;
 491 }