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 "AliExternalTrackParam.h"
  23 //#include "TMath.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
  43
  44   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
  45   Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
  46
  47   AliKFParticleBase::Initialize( Param, Cov, Charge, mass );
  48 }
  49
  50
  51 AliKFParticle::AliKFParticle( const AliExternalTrackParam &track, Int_t PID )
  52 {
  53   // Constructor from ALICE track, PID hypothesis should be provided
  54
  55   track.GetXYZ(fP);
  56   track.GetPxPyPz(fP+3);
  57   fQ = (track.GetSigned1Pt() >0 ) ?1 :-1;
  58   track.GetCovarianceXYZPxPyPz( fC );
  59   Create(fP,fC,fQ,PID);
  60 }
  61
  62 AliKFParticle::AliKFParticle( const AliESDVertex &vertex )
  63 {
  64   // Constructor from ALICE vertex
  65
  66   vertex.GetXYZ( fP );
  67   vertex.GetCovMatrix( 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
  77 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
  78 {
  79   // Conversion to AliExternalTrackParam parameterization
  80
  81   Double_t cosA = p.GetPx(), sinA = p.GetPy();
  82   Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
  83   Double_t pti = 0;
  84   if( pt<1.e-4 ){
  85     cosA = 1;
  86     sinA = 0;
  87   } else {
  88     pti = 1./pt;
  89     cosA*=pti;
  90     sinA*=pti;
  91   }
  92   Alpha = TMath::ATan2(sinA,cosA);
  93   X   = p.GetX()*cosA + p.GetY()*sinA;
  94   P[0]= p.GetY()*cosA - p.GetX()*sinA;
  95   P[1]= p.GetZ();
  96   P[2]= 0;
  97   P[3]= p.GetPz()*pti;
  98   P[4]= p.GetQ()*pti;
  99 }
 100
 101
 102
 103 Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
 104 {
 105   //* Calculate distance from vertex [cm] in XY-plane
 106
 107   Double_t mP[8], mC[36];
 108   Transport( GetDStoPoint(vtx), mP, mC );
 109   Double_t d[2]={ vtx[0]-mP[0], vtx[1]-mP[1] };
 110   Double_t dist =  TMath::Sqrt( d[0]*d[0]+d[1]*d[1] );
 111   Double_t sign = d[0]*mP[3] - d[1]*mP[4];
 112   return (sign>=0) ?dist :-dist;
 113 }
 114
 115 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const
 116 {
 117   //* Calculate distance from vertex [cm] in XY-plane
 118
 119   return GetDistanceFromVertexXY( Vtx.fP );
 120 }
 121
 122 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliESDVertex &Vtx ) const
 123 {
 124   //* Calculate distance from vertex [cm] in XY-plane
 125
 126   return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
 127 }
 128
 129 Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const
 130 {
 131   //* Calculate distance to other particle [cm]
 132
 133   Double_t dS, dS1;
 134   GetDStoParticleXY( p, dS, dS1 );
 135   Double_t mP[8], mC[36], mP1[8], mC1[36];
 136   Transport( dS, mP, mC );
 137   p.Transport( dS1, mP1, mC1 );
 138   Double_t dx = mP[0]-mP1[0];
 139   Double_t dy = mP[1]-mP1[1];
 140   return TMath::Sqrt(dx*dx+dy*dy);
 141 }
 142
 143 Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const
 144 {
 145   //* Calculate sqrt(Chi2/ndf) deviation from other particle
 146
 147   Double_t dS, dS1;
 148   GetDStoParticleXY( p, dS, dS1 );
 149   Double_t mP1[8], mC1[36];
 150   p.Transport( dS1, mP1, mC1 );
 151
 152   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
 153
 154   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
 155                                         (mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );
 156
 157   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };
 158
 159   mC1[0] +=h[0]*h[0];
 160   mC1[1] +=h[1]*h[0];
 161   mC1[2] +=h[1]*h[1];
 162
 163   return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
 164 }
 165
 166
 167 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t v[], const Double_t Cv[] ) const
 168 {
 169   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 170   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 171
 172   Double_t mP[8];
 173   Double_t mC[36];
 174
 175   Transport( GetDStoPoint(v), mP, mC );
 176
 177   Double_t d[2]={ v[0]-mP[0], v[1]-mP[1] };
 178
 179   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
 180                                         (mP[3]*mP[3]+mP[4]*mP[4] )  );
 181
 182   Double_t h[2] = { mP[3]*sigmaS, mP[4]*sigmaS };
 183
 184   Double_t mSi[3] = { mC[0] +h[0]*h[0],
 185                       mC[1] +h[1]*h[0], mC[2] +h[1]*h[1] };
 186
 187   if( Cv ){
 188     mSi[0]+=Cv[0];
 189     mSi[1]+=Cv[1];
 190     mSi[2]+=Cv[2];
 191   }
 192   Double_t s = ( mSi[0]*mSi[2] - mSi[1]*mSi[1] );
 193   s = ( s > 1.E-20 )  ?1./s :0;
 194
 195   Double_t mS[3] = { mSi[2],
 196                      -mSi[1], mSi[0] };
 197
 198   return TMath::Sqrt( TMath::Abs(s*( ( mS[0]*d[0] + mS[1]*d[1] )*d[0]
 199                                      +(mS[1]*d[0] + mS[2]*d[1] )*d[1] ))/1);
 200 }
 201
 202
 203 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
 204 {
 205   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 206   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 207
 208   return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
 209 }
 210
 211 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliESDVertex &Vtx ) const
 212 {
 213   //* Calculate sqrt(Chi2/ndf) deviation from vertex
 214   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 215
 216   AliKFParticle v(Vtx);
 217   return GetDeviationFromVertexXY( v.fP, v.fC );
 218 }
 219
 220
 221 Double_t AliKFParticle::GetAngle  ( const AliKFParticle &p ) const
 222 {
 223   //* Calculate the opening angle between two particles
 224
 225   Double_t dS, dS1;
 226   GetDStoParticle( p, dS, dS1 );
 227   Double_t mP[8], mC[36], mP1[8], mC1[36];
 228   Transport( dS, mP, mC );
 229   p.Transport( dS1, mP1, mC1 );
 230   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
 231   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
 232   n*=n1;
 233   Double_t a = 0;
 234   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
 235   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 236   else a = (a>=0) ?0 :TMath::Pi();
 237   return a;
 238 }
 239
 240 Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const
 241 {
 242   //* Calculate the opening angle between two particles in XY plane
 243
 244   Double_t dS, dS1;
 245   GetDStoParticleXY( p, dS, dS1 );
 246   Double_t mP[8], mC[36], mP1[8], mC1[36];
 247   Transport( dS, mP, mC );
 248   p.Transport( dS1, mP1, mC1 );
 249   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 250   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
 251   n*=n1;
 252   Double_t a = 0;
 253   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
 254   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 255   else a = (a>=0) ?0 :TMath::Pi();
 256   return a;
 257 }
 258
 259 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const
 260 {
 261   //* Calculate the opening angle between two particles in RZ plane
 262
 263   Double_t dS, dS1;
 264   GetDStoParticle( p, dS, dS1 );
 265   Double_t mP[8], mC[36], mP1[8], mC1[36];
 266   Transport( dS, mP, mC );
 267   p.Transport( dS1, mP1, mC1 );
 268   Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
 269   Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
 270   Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
 271   Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
 272   n*=n1;
 273   Double_t a = 0;
 274   if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n;
 275   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
 276   else a = (a>=0) ?0 :TMath::Pi();
 277   return a;
 278 }