STEER/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 ClassImp(AliKFParticle)
  26
  27 Double_t AliKFParticle::fgBz = -5.;  //* Bz compoment of the magnetic field
  28
  29 void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
  30 {
  31   // Constructor from "cartesian" track, PID hypothesis should be provided
  32   //
  33   // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
  34   // Cov [21] = lower-triangular part of the covariance matrix:
  35   //
  36   //                (  0  .  .  .  .  . )
  37   //                (  1  2  .  .  .  . )
  38   //  Cov. matrix = (  3  4  5  .  .  . ) - numbering of covariance elements in Cov[]
  39   //                (  6  7  8  9  .  . )
  40   //                ( 10 11 12 13 14  . )
  41   //                ( 15 16 17 18 19 20 )
  42   Double_t C[21];
  43   for( int i=0; i<21; i++ ) C[i] = Cov[i];
  44
  45   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
  46   Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
  47
  48   AliKFParticleBase::Initialize( Param, C, Charge, mass );
  49 }
  50
  51 AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
  52 {
  53   // Constructor from ALICE track, PID hypothesis should be provided
  54
  55   track.GetXYZ(fP);
  56   track.PxPyPz(fP+3);
  57   fQ = track.Charge();
  58   track.GetCovarianceXYZPxPyPz( fC );
  59   Create(fP,fC,fQ,PID);
  60 }
  61
  62 AliKFParticle::AliKFParticle( const AliVVertex &vertex )
  63 {
  64   // Constructor from ALICE vertex
  65
  66   vertex.GetXYZ( fP );
  67   vertex.GetCovarianceMatrix( fC );
  68   fChi2 = vertex.GetChi2();
  69   fNDF = 2*vertex.GetNContributors() - 3;
  70   fQ = 0;
  71   fAtProductionVertex = 0;
  72   fIsLinearized = 0;
  73   fSFromDecay = 0;
  74 }
  75
  76 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
  77 {
  78   // Conversion to AliExternalTrackParam parameterization
  79
  80   Double_t cosA = p.GetPx(), sinA = p.GetPy();
  81   Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
  82   Double_t pti = 0;
  83   if( pt<1.e-4 ){
  84     cosA = 1;
  85     sinA = 0;
  86   } else {
  87     pti = 1./pt;
  88     cosA*=pti;
  89     sinA*=pti;
  90   }
  91   Alpha = TMath::ATan2(sinA,cosA);
  92   X   = p.GetX()*cosA + p.GetY()*sinA;
  93   P[0]= p.GetY()*cosA - p.GetX()*sinA;
  94   P[1]= p.GetZ();
  95   P[2]= 0;
  96   P[3]= p.GetPz()*pti;
  97   P[4]= p.GetQ()*pti;
  98 }
  99
 100 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
 101 {
 102   //* Calculate DCA distance from vertex (transverse impact parameter) in XY
 103   //* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix
 104
 105   Bool_t ret = 0;
 106
 107   Double_t mP[8];
 108   Double_t mC[36];
 109
 110   Transport( GetDStoPoint(vtx), mP, mC );
 111
 112   Double_t dx = mP[0] - vtx[0];
 113   Double_t dy = mP[1] - vtx[1];
 114   Double_t px = mP[3];
 115   Double_t py = mP[4];
 116   Double_t pt = TMath::Sqrt(px*px + py*py);
 117   Double_t ex=0, ey=0;
 118   if( pt<1.e-4 ){
 119     ret = 1;
 120     pt = 1.;
 121     val = 1.e4;
 122   } else{
 123     ex = px/pt;
 124     ey = py/pt;
 125     val = dy*ex - dx*ey;
 126   }
 127
 128   Double_t h0 = -ey;
 129   Double_t h1 = ex;
 130   Double_t h3 = (dy*ey + dx*ex)*ey/pt;
 131   Double_t h4 = -(dy*ey + dx*ex)*ex/pt;
 132
 133   err =
 134     h0*(h0*GetCovariance(0,0) + h1*GetCovariance(0,1) + h3*GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
 135     h1*(h0*GetCovariance(1,0) + h1*GetCovariance(1,1) + h3*GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
 136     h3*(h0*GetCovariance(3,0) + h1*GetCovariance(3,1) + h3*GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
 137     h4*(h0*GetCovariance(4,0) + h1*GetCovariance(4,1) + h3*GetCovariance(4,3) + h4*GetCovariance(4,4) );
 138
 139   if( Cv ){
 140     err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] );
 141   }
 142
 143   err = TMath::Sqrt(TMath::Abs(err));
 144
 145   return ret;
 146 }
 147
 148 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
 149 {
 150   return GetDistanceFromVertexXY( vtx, 0, val, err );
 151 }
 152
 153
 154 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const
 155 {
 156   //* Calculate distance from vertex [cm] in XY-plane
 157
 158   return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
 159 }
 160
 161 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx, Double_t &val, Double_t &err ) const
 162 {
 163   //* Calculate distance from vertex [cm] in XY-plane
 164
 165   return GetDistanceFromVertexXY( AliKFParticle(Vtx), val, err );
 166 }
 167
 168 Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
 169 {
 170   //* Calculate distance from vertex [cm] in XY-plane
 171   Double_t val, err;
 172   GetDistanceFromVertexXY( vtx, 0, val, err );
 173   return val;
 174 }
 175
 176 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const
 177 {
 178   //* Calculate distance from vertex [cm] in XY-plane
 179
 180   return GetDistanceFromVertexXY( Vtx.fP );
 181 }
 182
 183 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const
 184 {
 185   //* Calculate distance from vertex [cm] in XY-plane
 186
 187   return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
 188 }
 189
 190 Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const
 191 {
 192   //* Calculate distance to other particle [cm]
 193
 194   Double_t dS, dS1;
 195   GetDStoParticleXY( p, dS, dS1 );
 196   Double_t mP[8], mC[36], mP1[8], mC1[36];
 197   Transport( dS, mP, mC );
 198   p.Transport( dS1, mP1, mC1 );
 199   Double_t dx = mP[0]-mP1[0];
 200   Double_t dy = mP[1]-mP1[1];
 201   return TMath::Sqrt(dx*dx+dy*dy);
 202 }
 203
 204 Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const
 205 {
 206   //* Calculate sqrt(Chi2/ndf) deviation from other particle
 207
 208   Double_t dS, dS1;
 209   GetDStoParticleXY( p, dS, dS1 );
 210   Double_t mP1[8], mC1[36];
 211   p.Transport( dS1, mP1, mC1 );
 212
 213   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
 214
 215   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
 216                                         (mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );
 217
 218   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };
 219
 220   mC1[0] +=h[0]*h[0];
 221   mC1[1] +=h[1]*h[0];
 222   mC1[2] +=h[1]*h[1];
 223
 224   return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
 225 }
 226
 227
 228 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const
 229 {
 230   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 231   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 232
 233   Double_t val, err;
 234   Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
 235   if( problem || err<1.e-20 ) return 1.e4;
 236   else return val/err;
 237 }
 238
 239
 240 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
 241 {
 242   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 243   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 244
 245   return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
 246 }
 247
 248 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const
 249 {
 250   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 251   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 252
 253   AliKFParticle v(Vtx);
 254   return GetDeviationFromVertexXY( v.fP, v.fC );
 255 }
 256
 257 Double_t AliKFParticle::GetAngle  ( const AliKFParticle &p ) const
 258 {
 259   //* Calculate the opening angle between two particles
 260
 261   Double_t dS, dS1;
 262   GetDStoParticle( p, dS, dS1 );
 263   Double_t mP[8], mC[36], mP1[8], mC1[36];
 264   Transport( dS, mP, mC );
 265   p.Transport( dS1, mP1, mC1 );
 266   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
 267   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
 268   n*=n1;
 269   Double_t a = 0;
 270   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
 271   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 272   else a = (a>=0) ?0 :TMath::Pi();
 273   return a;
 274 }
 275
 276 Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const
 277 {
 278   //* Calculate the opening angle between two particles in XY plane
 279
 280   Double_t dS, dS1;
 281   GetDStoParticleXY( p, dS, dS1 );
 282   Double_t mP[8], mC[36], mP1[8], mC1[36];
 283   Transport( dS, mP, mC );
 284   p.Transport( dS1, mP1, mC1 );
 285   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 286   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
 287   n*=n1;
 288   Double_t a = 0;
 289   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
 290   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 291   else a = (a>=0) ?0 :TMath::Pi();
 292   return a;
 293 }
 294
 295 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const
 296 {
 297   //* Calculate the opening angle between two particles in RZ plane
 298
 299   Double_t dS, dS1;
 300   GetDStoParticle( p, dS, dS1 );
 301   Double_t mP[8], mC[36], mP1[8], mC1[36];
 302   Transport( dS, mP, mC );
 303   p.Transport( dS1, mP1, mC1 );
 304   Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 305   Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
 306   Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
 307   Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
 308   n*=n1;
 309   Double_t a = 0;
 310   if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n;
 311   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 312   else a = (a>=0) ?0 :TMath::Pi();
 313   return a;
 314 }
 315
 316
 317 /*
 318
 319 #include "AliExternalTrackParam.h"
 320
 321 void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p,
 322                                            Double_t &DS, Double_t &DS1 )
 323   const
 324 {
 325   DS = DS1 = 0;
 326   Double_t x1, a1, x2, a2;
 327   Double_t par1[5], par2[5], cov[15];
 328   for(int i=0; i<15; i++) cov[i] = 0;
 329   cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;
 330
 331   GetExternalTrackParam( *this, x1, a1, par1 );
 332   GetExternalTrackParam( p, x2, a2, par2 );
 333
 334   AliExternalTrackParam t1(x1,a1, par1, cov);
 335   AliExternalTrackParam t2(x2,a2, par2, cov);
 336
 337   Double_t xe1=0, xe2=0;
 338   t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
 339   t1.PropagateTo( xe1, -GetFieldAlice() );
 340   t2.PropagateTo( xe2, -GetFieldAlice() );
 341
 342   Double_t xyz1[3], xyz2[3];
 343   t1.GetXYZ( xyz1 );
 344   t2.GetXYZ( xyz2 );
 345
 346   DS = GetDStoPoint( xyz1 );
 347   DS1 = p.GetDStoPoint( xyz2 );
 348
 349   return;
 350 }
 351 */