1 //----------------------------------------------------------------------------
2 // Implementation of the AliKFParticle class
4 // @author S.Gorbunov, I.Kisel
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
13 // This class is ALICE interface to general mathematics in AliKFParticleCore
15 // -= Copyright © ALICE HLT Group =-
16 //____________________________________________________________________________
19 #include "AliKFParticle.h"
20 #include "TDatabasePDG.h"
21 #include "TParticlePDG.h"
22 #include "AliVTrack.h"
23 #include "AliVVertex.h"
25 #include "AliExternalTrackParam.h"
27 ClassImp(AliKFParticle)
29 Double_t AliKFParticle::fgBz = -5.; //* Bz compoment of the magnetic field
31 AliKFParticle::AliKFParticle( const AliKFParticle &d1, const AliKFParticle &d2, Bool_t gamma )
39 ConstructGamma(d1, d2);
42 void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
44 // Constructor from "cartesian" track, PID hypothesis should be provided
46 // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
47 // Cov [21] = lower-triangular part of the covariance matrix:
51 // Cov. matrix = ( 3 4 5 . . . ) - numbering of covariance elements in Cov[]
53 // ( 10 11 12 13 14 . )
54 // ( 15 16 17 18 19 20 )
56 for( int i=0; i<21; i++ ) C[i] = Cov[i];
58 TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
59 Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
61 AliKFParticleBase::Initialize( Param, C, Charge, mass );
64 void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Double_t Mass )
66 // Constructor from "cartesian" track, PID hypothesis should be provided
68 // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
69 // Cov [21] = lower-triangular part of the covariance matrix:
73 // Cov. matrix = ( 3 4 5 . . . ) - numbering of covariance elements in Cov[]
75 // ( 10 11 12 13 14 . )
76 // ( 15 16 17 18 19 20 )
78 for( int i=0; i<21; i++ ) C[i] = Cov[i];
80 AliKFParticleBase::Initialize( Param, C, Charge, Mass );
83 AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
85 // Constructor from ALICE track, PID hypothesis should be provided
90 track.GetCovarianceXYZPxPyPz( fC );
94 AliKFParticle::AliKFParticle( const AliExternalTrackParam &track, Double_t Mass, Int_t Charge )
96 // Constructor from ALICE track, Mass and Charge hypothesis should be provided
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);
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]));
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];
114 m43*=track.GetSign();
115 m44*=track.GetSign();
116 m45*=track.GetSign();
118 const Double_t *cTr = track.GetCovariance();
120 fC[0 ] = cTr[0]*m00*m00;
121 fC[1 ] = cTr[0]*m00*m10;
122 fC[2 ] = cTr[0]*m10*m10;
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);
142 Create(fP,fC,fQ,Mass);
145 AliKFParticle::AliKFParticle( const AliVVertex &vertex )
147 // Constructor from ALICE vertex
150 vertex.GetCovarianceMatrix( fC );
151 fChi2 = vertex.GetChi2();
152 fNDF = 2*vertex.GetNContributors() - 3;
154 fAtProductionVertex = 0;
159 void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
161 // Conversion to AliExternalTrackParam parameterization
163 Double_t cosA = p.GetPx(), sinA = p.GetPy();
164 Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
174 Alpha = TMath::ATan2(sinA,cosA);
175 X = p.GetX()*cosA + p.GetY()*sinA;
176 P[0]= p.GetY()*cosA - p.GetX()*sinA;
183 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
185 //* Calculate DCA distance from vertex (transverse impact parameter) in XY
186 //* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix
193 Transport( GetDStoPoint(vtx), mP, mC );
195 Double_t dx = mP[0] - vtx[0];
196 Double_t dy = mP[1] - vtx[1];
199 Double_t pt = TMath::Sqrt(px*px + py*py);
213 Double_t h3 = (dy*ey + dx*ex)*ey/pt;
214 Double_t h4 = -(dy*ey + dx*ex)*ex/pt;
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) );
223 err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] );
226 err = TMath::Sqrt(TMath::Abs(err));
231 Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
233 return GetDistanceFromVertexXY( vtx, 0, val, err );
237 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const
239 //* Calculate distance from vertex [cm] in XY-plane
241 return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
244 Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx, Double_t &val, Double_t &err ) const
246 //* Calculate distance from vertex [cm] in XY-plane
248 return GetDistanceFromVertexXY( AliKFParticle(Vtx), val, err );
251 Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
253 //* Calculate distance from vertex [cm] in XY-plane
255 GetDistanceFromVertexXY( vtx, 0, val, err );
259 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const
261 //* Calculate distance from vertex [cm] in XY-plane
263 return GetDistanceFromVertexXY( Vtx.fP );
266 Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const
268 //* Calculate distance from vertex [cm] in XY-plane
270 return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
273 Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const
275 //* Calculate distance to other particle [cm]
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);
287 Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const
289 //* Calculate sqrt(Chi2/ndf) deviation from other particle
292 GetDStoParticleXY( p, dS, dS1 );
293 Double_t mP1[8], mC1[36];
294 p.Transport( dS1, mP1, mC1 );
296 Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
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] ) );
301 Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };
307 return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
311 Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const
313 //* Calculate sqrt(Chi2/ndf) deviation from vertex
314 //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
317 Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
318 if( problem || err<1.e-20 ) return 1.e4;
323 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
325 //* Calculate sqrt(Chi2/ndf) deviation from vertex
326 //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
328 return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
331 Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const
333 //* Calculate sqrt(Chi2/ndf) deviation from vertex
334 //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
336 AliKFParticle v(Vtx);
337 return GetDeviationFromVertexXY( v.fP, v.fC );
340 Double_t AliKFParticle::GetAngle ( const AliKFParticle &p ) const
342 //* Calculate the opening angle between two particles
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] );
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();
359 Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const
361 //* Calculate the opening angle between two particles in XY plane
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] );
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();
378 Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const
380 //* Calculate the opening angle between two particles in RZ plane
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] );
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();
402 #include "AliExternalTrackParam.h"
404 void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p,
405 Double_t &DS, Double_t &DS1 )
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;
414 GetExternalTrackParam( *this, x1, a1, par1 );
415 GetExternalTrackParam( p, x2, a2, par2 );
417 AliExternalTrackParam t1(x1,a1, par1, cov);
418 AliExternalTrackParam t2(x2,a2, par2, cov);
420 Double_t xe1=0, xe2=0;
421 t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
422 t1.PropagateTo( xe1, -GetFieldAlice() );
423 t2.PropagateTo( xe2, -GetFieldAlice() );
425 Double_t xyz1[3], xyz2[3];
429 DS = GetDStoPoint( xyz1 );
430 DS1 = p.GetDStoPoint( xyz2 );
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();
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,
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),
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;
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);
490 return ( dx*Px() + dy*Py() )*mipt2;