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