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
  52 AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
  53 {
  54   // Constructor from ALICE track, PID hypothesis should be provided
  55
  56   track.XvYvZv(fP);
  57   track.PxPyPz(fP+3);
  58   fQ = track.Charge();
  59   track.GetCovarianceXYZPxPyPz( fC );
  60   Create(fP,fC,fQ,PID);
  61 }
  62
  63 AliKFParticle::AliKFParticle( const AliVVertex &vertex )
  64 {
  65   // Constructor from ALICE vertex
  66
  67   vertex.GetXYZ( fP );
  68   vertex.GetCovarianceMatrix( fC );
  69   fChi2 = vertex.GetChi2();
  70   fNDF = 2*vertex.GetNContributors() - 3;
  71   fQ = 0;
  72   fAtProductionVertex = 0;
  73   fIsLinearized = 0;
  74   fSFromDecay = 0;
  75 }
  76
  77
  78 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
  79 {
  80   // Conversion to AliExternalTrackParam parameterization
  81
  82   Double_t cosA = p.GetPx(), sinA = p.GetPy();
  83   Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
  84   Double_t pti = 0;
  85   if( pt<1.e-4 ){
  86     cosA = 1;
  87     sinA = 0;
  88   } else {
  89     pti = 1./pt;
  90     cosA*=pti;
  91     sinA*=pti;
  92   }
  93   Alpha = TMath::ATan2(sinA,cosA);
  94   X   = p.GetX()*cosA + p.GetY()*sinA;
  95   P[0]= p.GetY()*cosA - p.GetX()*sinA;
  96   P[1]= p.GetZ();
  97   P[2]= 0;
  98   P[3]= p.GetPz()*pti;
  99   P[4]= p.GetQ()*pti;
 100 }
 101
 102
 103
 104 Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
 105 {
 106   //* Calculate distance from vertex [cm] in XY-plane
 107
 108   Double_t mP[8], mC[36];
 109   Transport( GetDStoPoint(vtx), mP, mC );
 110   Double_t d[2]={ vtx[0]-mP[0], vtx[1]-mP[1] };
 111   Double_t dist =  TMath::Sqrt( d[0]*d[0]+d[1]*d[1] );
 112   Double_t sign = d[0]*mP[3] - d[1]*mP[4];
 113   return (sign>=0) ?dist :-dist;
 114 }
 115
 116 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const
 117 {
 118   //* Calculate distance from vertex [cm] in XY-plane
 119
 120   return GetDistanceFromVertexXY( Vtx.fP );
 121 }
 122
 123 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const
 124 {
 125   //* Calculate distance from vertex [cm] in XY-plane
 126
 127   return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
 128 }
 129
 130 Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const
 131 {
 132   //* Calculate distance to other particle [cm]
 133
 134   Double_t dS, dS1;
 135   GetDStoParticleXY( p, dS, dS1 );
 136   Double_t mP[8], mC[36], mP1[8], mC1[36];
 137   Transport( dS, mP, mC );
 138   p.Transport( dS1, mP1, mC1 );
 139   Double_t dx = mP[0]-mP1[0];
 140   Double_t dy = mP[1]-mP1[1];
 141   return TMath::Sqrt(dx*dx+dy*dy);
 142 }
 143
 144 Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const
 145 {
 146   //* Calculate sqrt(Chi2/ndf) deviation from other particle
 147
 148   Double_t dS, dS1;
 149   GetDStoParticleXY( p, dS, dS1 );
 150   Double_t mP1[8], mC1[36];
 151   p.Transport( dS1, mP1, mC1 );
 152
 153   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
 154
 155   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
 156                                         (mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );
 157
 158   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };
 159
 160   mC1[0] +=h[0]*h[0];
 161   mC1[1] +=h[1]*h[0];
 162   mC1[2] +=h[1]*h[1];
 163
 164   return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
 165 }
 166
 167
 168 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t v[], const Double_t Cv[] ) const
 169 {
 170   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 171   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 172
 173   Double_t mP[8];
 174   Double_t mC[36];
 175
 176   Transport( GetDStoPoint(v), mP, mC );
 177
 178   Double_t d[2]={ v[0]-mP[0], v[1]-mP[1] };
 179
 180   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
 181                                         (mP[3]*mP[3]+mP[4]*mP[4] )  );
 182
 183   Double_t h[2] = { mP[3]*sigmaS, mP[4]*sigmaS };
 184
 185   Double_t mSi[3] = { mC[0] +h[0]*h[0],
 186                       mC[1] +h[1]*h[0], mC[2] +h[1]*h[1] };
 187
 188   if( Cv ){
 189     mSi[0]+=Cv[0];
 190     mSi[1]+=Cv[1];
 191     mSi[2]+=Cv[2];
 192   }
 193   Double_t s = ( mSi[0]*mSi[2] - mSi[1]*mSi[1] );
 194   s = ( s > 1.E-20 )  ?1./s :0;
 195
 196   Double_t mS[3] = { mSi[2],
 197                      -mSi[1], mSi[0] };
 198
 199   return TMath::Sqrt( TMath::Abs(s*( ( mS[0]*d[0] + mS[1]*d[1] )*d[0]
 200                                      +(mS[1]*d[0] + mS[2]*d[1] )*d[1] ))/1);
 201 }
 202
 203
 204 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
 205 {
 206   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 207   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 208
 209   return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
 210 }
 211
 212 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const
 213 {
 214   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 215   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 216
 217   AliKFParticle v(Vtx);
 218   return GetDeviationFromVertexXY( v.fP, v.fC );
 219 }
 220
 221
 222 Double_t AliKFParticle::GetAngle  ( const AliKFParticle &p ) const
 223 {
 224   //* Calculate the opening angle between two particles
 225
 226   Double_t dS, dS1;
 227   GetDStoParticle( p, dS, dS1 );
 228   Double_t mP[8], mC[36], mP1[8], mC1[36];
 229   Transport( dS, mP, mC );
 230   p.Transport( dS1, mP1, mC1 );
 231   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
 232   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
 233   n*=n1;
 234   Double_t a = 0;
 235   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
 236   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 237   else a = (a>=0) ?0 :TMath::Pi();
 238   return a;
 239 }
 240
 241 Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const
 242 {
 243   //* Calculate the opening angle between two particles in XY plane
 244
 245   Double_t dS, dS1;
 246   GetDStoParticleXY( p, dS, dS1 );
 247   Double_t mP[8], mC[36], mP1[8], mC1[36];
 248   Transport( dS, mP, mC );
 249   p.Transport( dS1, mP1, mC1 );
 250   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 251   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
 252   n*=n1;
 253   Double_t a = 0;
 254   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
 255   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 256   else a = (a>=0) ?0 :TMath::Pi();
 257   return a;
 258 }
 259
 260 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const
 261 {
 262   //* Calculate the opening angle between two particles in RZ plane
 263
 264   Double_t dS, dS1;
 265   GetDStoParticle( p, dS, dS1 );
 266   Double_t mP[8], mC[36], mP1[8], mC1[36];
 267   Transport( dS, mP, mC );
 268   p.Transport( dS1, mP1, mC1 );
 269   Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 270   Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
 271   Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
 272   Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
 273   n*=n1;
 274   Double_t a = 0;
 275   if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n;
 276   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 277   else a = (a>=0) ?0 :TMath::Pi();
 278   return a;
 279 }
 280
 281
 282 /*
 283
 284 #include "AliExternalTrackParam.h"
 285
 286 void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p,
 287                                            Double_t &DS, Double_t &DS1 )
 288   const
 289 {
 290   DS = DS1 = 0;
 291   Double_t x1, a1, x2, a2;
 292   Double_t par1[5], par2[5], cov[15];
 293   for(int i=0; i<15; i++) cov[i] = 0;
 294   cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;
 295
 296   GetExternalTrackParam( *this, x1, a1, par1 );
 297   GetExternalTrackParam( p, x2, a2, par2 );
 298
 299   AliExternalTrackParam t1(x1,a1, par1, cov);
 300   AliExternalTrackParam t2(x2,a2, par2, cov);
 301
 302   Double_t xe1=0, xe2=0;
 303   t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
 304   t1.PropagateTo( xe1, -GetFieldAlice() );
 305   t2.PropagateTo( xe2, -GetFieldAlice() );
 306
 307   Double_t xyz1[3], xyz2[3];
 308   t1.GetXYZ( xyz1 );
 309   t2.GetXYZ( xyz2 );
 310
 311   DS = GetDStoPoint( xyz1 );
 312   DS1 = p.GetDStoPoint( xyz2 );
 313
 314   return;
 315 }
 316 */