[u/mrichter/AliRoot.git] / STEER / ESD / AliKFParticle.cxx

//----------------------------------------------------------------------------
// Implementation of the AliKFParticle class
// .
// @author  S.Gorbunov, I.Kisel
// @version 1.0
// @since   13.05.07
// 
// Class to reconstruct and store the decayed particle parameters.
// The method is described in CBM-SOFT note 2007-003, 
// ``Reconstruction of decayed particles based on the Kalman filter'', 
// http://www.gsi.de/documents/DOC-2007-May-14-1.pdf
//
// This class is ALICE interface to general mathematics in AliKFParticleCore
// 
//  -= Copyright &copy ALICE HLT Group =-
//____________________________________________________________________________


#include "AliKFParticle.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"
#include "AliVTrack.h"
#include "AliVVertex.h"

#include "AliExternalTrackParam.h"

ClassImp(AliKFParticle)

Double_t AliKFParticle::fgBz = -5.;  //* Bz compoment of the magnetic field

AliKFParticle::AliKFParticle( const AliKFParticle &d1, const AliKFParticle &d2, Bool_t gamma )
{
  if (!gamma) {
    AliKFParticle mother;
    mother+= d1;
    mother+= d2;
    *this = mother;
  } else
    ConstructGamma(d1, d2);
}

void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
{
  // Constructor from "cartesian" track, PID hypothesis should be provided
  //
  // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
  // Cov [21] = lower-triangular part of the covariance matrix:
  //
  //                (  0  .  .  .  .  . )
  //                (  1  2  .  .  .  . )
  //  Cov. matrix = (  3  4  5  .  .  . ) - numbering of covariance elements in Cov[]
  //                (  6  7  8  9  .  . )
  //                ( 10 11 12 13 14  . )
  //                ( 15 16 17 18 19 20 )
  Double_t C[21];
  for( int i=0; i<21; i++ ) C[i] = Cov[i];
  
  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
  Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
  
  AliKFParticleBase::Initialize( Param, C, Charge, mass );
}

void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Double_t Mass )
{
  // Constructor from "cartesian" track, PID hypothesis should be provided
  //
  // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
  // Cov [21] = lower-triangular part of the covariance matrix:
  //
  //                (  0  .  .  .  .  . )
  //                (  1  2  .  .  .  . )
  //  Cov. matrix = (  3  4  5  .  .  . ) - numbering of covariance elements in Cov[]
  //                (  6  7  8  9  .  . )
  //                ( 10 11 12 13 14  . )
  //                ( 15 16 17 18 19 20 )
  Double_t C[21];
  for( int i=0; i<21; i++ ) C[i] = Cov[i];

  AliKFParticleBase::Initialize( Param, C, Charge, Mass );
}

AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
{
  // Constructor from ALICE track, PID hypothesis should be provided

  track.GetXYZ(fP);
  track.PxPyPz(fP+3);
  fQ = track.Charge();
  track.GetCovarianceXYZPxPyPz( fC );
  Create(fP,fC,fQ,PID);
}

AliKFParticle::AliKFParticle( const AliExternalTrackParam &track, Double_t Mass, Int_t Charge )
{
  // Constructor from ALICE track, Mass and Charge hypothesis should be provided

  track.GetXYZ(fP);
  track.PxPyPz(fP+3);
  fQ = track.Charge()*TMath::Abs(Charge);
  fP[3] *= TMath::Abs(Charge);
  fP[4] *= TMath::Abs(Charge);
  fP[5] *= TMath::Abs(Charge);

  Double_t pt=1./TMath::Abs(track.GetParameter()[4]) * TMath::Abs(Charge);
  Double_t cs=TMath::Cos(track.GetAlpha()), sn=TMath::Sin(track.GetAlpha());
  Double_t r=TMath::Sqrt((1.-track.GetParameter()[2])*(1.+track.GetParameter()[2]));

  Double_t m00=-sn, m10=cs;
  Double_t m23=-pt*(sn + track.GetParameter()[2]*cs/r), m43=-pt*pt*(r*cs - track.GetParameter()[2]*sn);
  Double_t m24= pt*(cs - track.GetParameter()[2]*sn/r), m44=-pt*pt*(r*sn + track.GetParameter()[2]*cs);
  Double_t m35=pt, m45=-pt*pt*track.GetParameter()[3];

  m43*=track.GetSign();
  m44*=track.GetSign();
  m45*=track.GetSign();

  const Double_t *cTr = track.GetCovariance();

  fC[0 ] = cTr[0]*m00*m00;
  fC[1 ] = cTr[0]*m00*m10; 
  fC[2 ] = cTr[0]*m10*m10;
  fC[3 ] = cTr[1]*m00; 
  fC[4 ] = cTr[1]*m10; 
  fC[5 ] = cTr[2];
  fC[6 ] = m00*(cTr[3]*m23 + cTr[10]*m43); 
  fC[7 ] = m10*(cTr[3]*m23 + cTr[10]*m43); 
  fC[8 ] = cTr[4]*m23 + cTr[11]*m43; 
  fC[9 ] = m23*(cTr[5]*m23 + cTr[12]*m43)  +  m43*(cTr[12]*m23 + cTr[14]*m43);
  fC[10] = m00*(cTr[3]*m24 + cTr[10]*m44); 
  fC[11] = m10*(cTr[3]*m24 + cTr[10]*m44); 
  fC[12] = cTr[4]*m24 + cTr[11]*m44; 
  fC[13] = m23*(cTr[5]*m24 + cTr[12]*m44)  +  m43*(cTr[12]*m24 + cTr[14]*m44);
  fC[14] = m24*(cTr[5]*m24 + cTr[12]*m44)  +  m44*(cTr[12]*m24 + cTr[14]*m44);
  fC[15] = m00*(cTr[6]*m35 + cTr[10]*m45); 
  fC[16] = m10*(cTr[6]*m35 + cTr[10]*m45); 
  fC[17] = cTr[7]*m35 + cTr[11]*m45; 
  fC[18] = m23*(cTr[8]*m35 + cTr[12]*m45)  +  m43*(cTr[13]*m35 + cTr[14]*m45);
  fC[19] = m24*(cTr[8]*m35 + cTr[12]*m45)  +  m44*(cTr[13]*m35 + cTr[14]*m45); 
  fC[20] = m35*(cTr[9]*m35 + cTr[13]*m45)  +  m45*(cTr[13]*m35 + cTr[14]*m45);

  Create(fP,fC,fQ,Mass);
}

AliKFParticle::AliKFParticle( const AliVVertex &vertex )
{
  // Constructor from ALICE vertex

  vertex.GetXYZ( fP );
  vertex.GetCovarianceMatrix( fC );  
  fChi2 = vertex.GetChi2();
  fNDF = 2*vertex.GetNContributors() - 3;
  fQ = 0;
  fAtProductionVertex = 0;
  fIsLinearized = 0;
  fSFromDecay = 0;
}

void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] ) 
{
  // Conversion to AliExternalTrackParam parameterization

  Double_t cosA = p.GetPx(), sinA = p.GetPy(); 
  Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
  Double_t pti = 0;
  if( pt<1.e-4 ){
    cosA = 1;
    sinA = 0;
  } else {
    pti = 1./pt;
    cosA*=pti;
    sinA*=pti;
  }
  Alpha = TMath::ATan2(sinA,cosA);  
  X   = p.GetX()*cosA + p.GetY()*sinA;   
  P[0]= p.GetY()*cosA - p.GetX()*sinA;
  P[1]= p.GetZ();
  P[2]= 0;
  P[3]= p.GetPz()*pti;
  P[4]= p.GetQ()*pti;
}

Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
{
  //* Calculate DCA distance from vertex (transverse impact parameter) in XY
  //* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix

  Bool_t ret = 0;
  
  Double_t mP[8];
  Double_t mC[36];
  
  Transport( GetDStoPoint(vtx), mP, mC );  

  Double_t dx = mP[0] - vtx[0];
  Double_t dy = mP[1] - vtx[1];
  Double_t px = mP[3];
  Double_t py = mP[4];
  Double_t pt = TMath::Sqrt(px*px + py*py);
  Double_t ex=0, ey=0;
  if( pt<1.e-4 ){
    ret = 1;
    pt = 1.;
    val = 1.e4;
  } else{
    ex = px/pt;
    ey = py/pt;
    val = dy*ex - dx*ey;
  }

  Double_t h0 = -ey;
  Double_t h1 = ex;
  Double_t h3 = (dy*ey + dx*ex)*ey/pt;
  Double_t h4 = -(dy*ey + dx*ex)*ex/pt;
  
  err = 
    h0*(h0*GetCovariance(0,0) + h1*GetCovariance(0,1) + h3*GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
    h1*(h0*GetCovariance(1,0) + h1*GetCovariance(1,1) + h3*GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
    h3*(h0*GetCovariance(3,0) + h1*GetCovariance(3,1) + h3*GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
    h4*(h0*GetCovariance(4,0) + h1*GetCovariance(4,1) + h3*GetCovariance(4,3) + h4*GetCovariance(4,4) );

  if( Cv ){
    err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] ); 
  }

  err = TMath::Sqrt(TMath::Abs(err));

  return ret;
}

Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
{
  return GetDistanceFromVertexXY( vtx, 0, val, err );
}


Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const 
{
  //* Calculate distance from vertex [cm] in XY-plane

  return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
}

Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx, Double_t &val, Double_t &err ) const 
{
  //* Calculate distance from vertex [cm] in XY-plane

  return GetDistanceFromVertexXY( AliKFParticle(Vtx), val, err );
}

Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
{
  //* Calculate distance from vertex [cm] in XY-plane
  Double_t val, err;
  GetDistanceFromVertexXY( vtx, 0, val, err );
  return val;
}

Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const 
{
  //* Calculate distance from vertex [cm] in XY-plane

  return GetDistanceFromVertexXY( Vtx.fP );
}

Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const 
{
  //* Calculate distance from vertex [cm] in XY-plane

  return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
}

Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const 
{
  //* Calculate distance to other particle [cm]

  Double_t dS, dS1;
  GetDStoParticleXY( p, dS, dS1 );   
  Double_t mP[8], mC[36], mP1[8], mC1[36];
  Transport( dS, mP, mC ); 
  p.Transport( dS1, mP1, mC1 ); 
  Double_t dx = mP[0]-mP1[0]; 
  Double_t dy = mP[1]-mP1[1]; 
  return TMath::Sqrt(dx*dx+dy*dy);
}

Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const 
{
  //* Calculate sqrt(Chi2/ndf) deviation from other particle

  Double_t dS, dS1;
  GetDStoParticleXY( p, dS, dS1 );   
  Double_t mP1[8], mC1[36];
  p.Transport( dS1, mP1, mC1 ); 

  Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };

  Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
					(mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );

  Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };       
  
  mC1[0] +=h[0]*h[0];
  mC1[1] +=h[1]*h[0]; 
  mC1[2] +=h[1]*h[1]; 

  return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
}


Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const 
{
  //* Calculate sqrt(Chi2/ndf) deviation from vertex
  //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix

  Double_t val, err;
  Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
  if( problem || err<1.e-20 ) return 1.e4;
  else return val/err;
}


Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const  
{
  //* Calculate sqrt(Chi2/ndf) deviation from vertex
  //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix

  return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
}

Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const 
{
  //* Calculate sqrt(Chi2/ndf) deviation from vertex
  //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix

  AliKFParticle v(Vtx);
  return GetDeviationFromVertexXY( v.fP, v.fC );
}

Double_t AliKFParticle::GetAngle  ( const AliKFParticle &p ) const 
{
  //* Calculate the opening angle between two particles

  Double_t dS, dS1;
  GetDStoParticle( p, dS, dS1 );   
  Double_t mP[8], mC[36], mP1[8], mC1[36];
  Transport( dS, mP, mC ); 
  p.Transport( dS1, mP1, mC1 ); 
  Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
  Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
  n*=n1;
  Double_t a = 0;
  if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
  if (TMath::Abs(a)<1.) a = TMath::ACos(a);
  else a = (a>=0) ?0 :TMath::Pi();
  return a;
}

Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const 
{
  //* Calculate the opening angle between two particles in XY plane

  Double_t dS, dS1;
  GetDStoParticleXY( p, dS, dS1 );   
  Double_t mP[8], mC[36], mP1[8], mC1[36];
  Transport( dS, mP, mC ); 
  p.Transport( dS1, mP1, mC1 ); 
  Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
  Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
  n*=n1;
  Double_t a = 0;
  if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
  if (TMath::Abs(a)<1.) a = TMath::ACos(a);
  else a = (a>=0) ?0 :TMath::Pi();
  return a;
}

Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const 
{
  //* Calculate the opening angle between two particles in RZ plane

  Double_t dS, dS1;
  GetDStoParticle( p, dS, dS1 );   
  Double_t mP[8], mC[36], mP1[8], mC1[36];
  Transport( dS, mP, mC ); 
  p.Transport( dS1, mP1, mC1 ); 
  Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
  Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
  Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
  Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
  n*=n1;
  Double_t a = 0;
  if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n; 
  if (TMath::Abs(a)<1.) a = TMath::ACos(a);
  else a = (a>=0) ?0 :TMath::Pi();
  return a;
}


/*

#include "AliExternalTrackParam.h"

void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p, 
					   Double_t &DS, Double_t &DS1 ) 
  const
{ 
  DS = DS1 = 0;   
  Double_t x1, a1, x2, a2;
  Double_t par1[5], par2[5], cov[15];
  for(int i=0; i<15; i++) cov[i] = 0;
  cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;

  GetExternalTrackParam( *this, x1, a1, par1 );
  GetExternalTrackParam( p, x2, a2, par2 );

  AliExternalTrackParam t1(x1,a1, par1, cov);
  AliExternalTrackParam t2(x2,a2, par2, cov);

  Double_t xe1=0, xe2=0;
  t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
  t1.PropagateTo( xe1, -GetFieldAlice() );
  t2.PropagateTo( xe2, -GetFieldAlice() );

  Double_t xyz1[3], xyz2[3];
  t1.GetXYZ( xyz1 );
  t2.GetXYZ( xyz2 );
  
  DS = GetDStoPoint( xyz1 );
  DS1 = p.GetDStoPoint( xyz2 );

  return;
}
*/

  // * Pseudo Proper Time of decay = (r*pt) / |pt| * M/|pt|
Double_t AliKFParticle::GetPseudoProperDecayTime( const AliKFParticle &pV, const Double_t& mass, Double_t* timeErr2 ) const
{ // TODO optimize with respect to time and stability
  const Double_t ipt2 = 1/( Px()*Px() + Py()*Py() );
  const Double_t mipt2 = mass*ipt2;
  const Double_t dx = X() - pV.X();
  const Double_t dy = Y() - pV.Y();

  if ( timeErr2 ) {
      // -- calculate error = sigma(f(r)) = f'Cf'
      // r = {x,y,px,py,x_pV,y_pV}
      // df/dr = { px*m/pt^2,
      //     py*m/pt^2,
      //    ( x - x_pV )*m*(1/pt^2 - 2(px/pt^2)^2),
      //    ( y - y_pV )*m*(1/pt^2 - 2(py/pt^2)^2),
      //     -px*m/pt^2,
      //     -py*m/pt^2 }
    const Double_t f0 = Px()*mipt2;
    const Double_t f1 = Py()*mipt2;
    const Double_t mipt2derivative = mipt2*(1-2*Px()*Px()*ipt2);
    const Double_t f2 = dx*mipt2derivative;
    const Double_t f3 = -dy*mipt2derivative;
    const Double_t f4 = -f0;
    const Double_t f5 = -f1;

    const Double_t& mC00 =    GetCovariance(0,0);
    const Double_t& mC10 =    GetCovariance(0,1);
    const Double_t& mC11 =    GetCovariance(1,1);
    const Double_t& mC20 =    GetCovariance(3,0);
    const Double_t& mC21 =    GetCovariance(3,1);
    const Double_t& mC22 =    GetCovariance(3,3);
    const Double_t& mC30 =    GetCovariance(4,0);
    const Double_t& mC31 =    GetCovariance(4,1);
    const Double_t& mC32 =    GetCovariance(4,3);
    const Double_t& mC33 =    GetCovariance(4,4);
    const Double_t& mC44 = pV.GetCovariance(0,0);
    const Double_t& mC54 = pV.GetCovariance(1,0);
    const Double_t& mC55 = pV.GetCovariance(1,1);

    *timeErr2 =
      f5*mC55*f5 +
      f5*mC54*f4 +
      f4*mC44*f4 +
      f3*mC33*f3 +
      f3*mC32*f2 +
      f3*mC31*f1 +
      f3*mC30*f0 +
      f2*mC22*f2 +
      f2*mC21*f1 +
      f2*mC20*f0 +
      f1*mC11*f1 +
      f1*mC10*f0 +
      f0*mC00*f0;
  }
  return ( dx*Px() + dy*Py() )*mipt2;
}
Commit	Line	Data
f826d409	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"
706952f5	22	#include "AliVTrack.h"
cfb5a2e1	23	#include "AliVVertex.h"
f826d409	24
706fa2e5	25	#include "AliExternalTrackParam.h"
706fa2e5	26
effa6338	27	ClassImp(AliKFParticle)
f826d409	28
616ffc76	29	Double_t AliKFParticle::fgBz = -5.; //* Bz compoment of the magnetic field
f826d409	30
537b06a4	31	AliKFParticle::AliKFParticle( const AliKFParticle &d1, const AliKFParticle &d2, Bool_t gamma )
	32	{
	33	if (!gamma) {
	34	AliKFParticle mother;
	35	mother+= d1;
	36	mother+= d2;
	37	*this = mother;
	38	} else
	39	ConstructGamma(d1, d2);
	40	}
	41
7972a8db	42	void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
f826d409	43	{
e7b09c95	44	// Constructor from "cartesian" track, PID hypothesis should be provided
	45	//
	46	// Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
	47	// Cov [21] = lower-triangular part of the covariance matrix:
	48	//
	49	// ( 0 . . . . . )
	50	// ( 1 2 . . . . )
	51	// Cov. matrix = ( 3 4 5 . . . ) - numbering of covariance elements in Cov[]
	52	// ( 6 7 8 9 . . )
	53	// ( 10 11 12 13 14 . )
	54	// ( 15 16 17 18 19 20 )
706952f5	55	Double_t C[21];
	56	for( int i=0; i<21; i++ ) C[i] = Cov[i];
	57
f826d409	58	TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
f826d409	59	Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
e7b09c95	60
706952f5	61	AliKFParticleBase::Initialize( Param, C, Charge, mass );
e7b09c95	62	}
616ffc76	63
73599358	64	void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Double_t Mass )
706fa2e5	65	{
	66	// Constructor from "cartesian" track, PID hypothesis should be provided
	67	//
	68	// Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
	69	// Cov [21] = lower-triangular part of the covariance matrix:
	70	//
	71	// ( 0 . . . . . )
	72	// ( 1 2 . . . . )
	73	// Cov. matrix = ( 3 4 5 . . . ) - numbering of covariance elements in Cov[]
	74	// ( 6 7 8 9 . . )
	75	// ( 10 11 12 13 14 . )
	76	// ( 15 16 17 18 19 20 )
	77	Double_t C[21];
	78	for( int i=0; i<21; i++ ) C[i] = Cov[i];
	79
	80	AliKFParticleBase::Initialize( Param, C, Charge, Mass );
	81	}
	82
706952f5	83	AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
e7b09c95	84	{
e7b09c95	85	// Constructor from ALICE track, PID hypothesis should be provided
616ffc76	86
b63c3605	87	track.GetXYZ(fP);
706952f5	88	track.PxPyPz(fP+3);
706952f5	89	fQ = track.Charge();
f826d409	90	track.GetCovarianceXYZPxPyPz( fC );
7972a8db	91	Create(fP,fC,fQ,PID);
f826d409	92	}
f826d409	93
73599358	94	AliKFParticle::AliKFParticle( const AliExternalTrackParam &track, Double_t Mass, Int_t Charge )
706fa2e5	95	{
	96	// Constructor from ALICE track, Mass and Charge hypothesis should be provided
	97
	98	track.GetXYZ(fP);
	99	track.PxPyPz(fP+3);
	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);
	104
	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]));
	108
	109	Double_t m00=-sn, m10=cs;
	110	Double_t m23=-pt(sn + track.GetParameter()[2]cs/r), m43=-ptpt(rcs - track.GetParameter()[2]sn);
	111	Double_t m24= pt(cs - track.GetParameter()[2]sn/r), m44=-ptpt(rsn + track.GetParameter()[2]cs);
	112	Double_t m35=pt, m45=-ptpttrack.GetParameter()[3];
	113
	114	m43*=track.GetSign();
	115	m44*=track.GetSign();
	116	m45*=track.GetSign();
	117
	118	const Double_t *cTr = track.GetCovariance();
	119
	120	fC[0 ] = cTr[0]m00m00;
	121	fC[1 ] = cTr[0]m00m10;
	122	fC[2 ] = cTr[0]m10m10;
	123	fC[3 ] = cTr[1]*m00;
	124	fC[4 ] = cTr[1]*m10;
	125	fC[5 ] = cTr[2];
	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);
	141
	142	Create(fP,fC,fQ,Mass);
	143	}
	144
706952f5	145	AliKFParticle::AliKFParticle( const AliVVertex &vertex )
f826d409	146	{
	147	// Constructor from ALICE vertex
	148
	149	vertex.GetXYZ( fP );
706952f5	150	vertex.GetCovarianceMatrix( fC );
f826d409	151	fChi2 = vertex.GetChi2();
	152	fNDF = 2*vertex.GetNContributors() - 3;
	153	fQ = 0;
	154	fAtProductionVertex = 0;
	155	fIsLinearized = 0;
	156	fSFromDecay = 0;
	157	}
616ffc76	158
616ffc76	159	void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
616ffc76	160	{
e7b09c95	161	// Conversion to AliExternalTrackParam parameterization
e7b09c95	162
616ffc76	163	Double_t cosA = p.GetPx(), sinA = p.GetPy();
	164	Double_t pt = TMath::Sqrt(cosAcosA + sinAsinA);
	165	Double_t pti = 0;
	166	if( pt<1.e-4 ){
	167	cosA = 1;
	168	sinA = 0;
	169	} else {
	170	pti = 1./pt;
	171	cosA*=pti;
	172	sinA*=pti;
	173	}
	174	Alpha = TMath::ATan2(sinA,cosA);
	175	X = p.GetX()cosA + p.GetY()sinA;
	176	P[0]= p.GetY()cosA - p.GetX()sinA;
	177	P[1]= p.GetZ();
	178	P[2]= 0;
	179	P[3]= p.GetPz()*pti;
	180	P[4]= p.GetQ()*pti;
	181	}
	182
446ce366	183	Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
	184	{
	185	//* Calculate DCA distance from vertex (transverse impact parameter) in XY
	186	//* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix
	187
	188	Bool_t ret = 0;
	189
	190	Double_t mP[8];
	191	Double_t mC[36];
	192
	193	Transport( GetDStoPoint(vtx), mP, mC );
	194
	195	Double_t dx = mP[0] - vtx[0];
	196	Double_t dy = mP[1] - vtx[1];
	197	Double_t px = mP[3];
	198	Double_t py = mP[4];
	199	Double_t pt = TMath::Sqrt(pxpx + pypy);
	200	Double_t ex=0, ey=0;
	201	if( pt<1.e-4 ){
	202	ret = 1;
	203	pt = 1.;
	204	val = 1.e4;
	205	} else{
	206	ex = px/pt;
924a4a7b	207	ey = py/pt;
446ce366	208	val = dyex - dxey;
446ce366	209	}
616ffc76	210
446ce366	211	Double_t h0 = -ey;
	212	Double_t h1 = ex;
	213	Double_t h3 = (dyey + dxex)*ey/pt;
	214	Double_t h4 = -(dyey + dxex)*ex/pt;
	215
	216	err =
	217	h0(h0GetCovariance(0,0) + h1GetCovariance(0,1) + h3GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
	218	h1(h0GetCovariance(1,0) + h1GetCovariance(1,1) + h3GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
	219	h3(h0GetCovariance(3,0) + h1GetCovariance(3,1) + h3GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
	220	h4(h0GetCovariance(4,0) + h1GetCovariance(4,1) + h3GetCovariance(4,3) + h4*GetCovariance(4,4) );
616ffc76	221
446ce366	222	if( Cv ){
	223	err+= h0(h0Cv[0] + h1Cv[1] ) + h1(h0Cv[1] + h1Cv[2] );
	224	}
	225
	226	err = TMath::Sqrt(TMath::Abs(err));
	227
	228	return ret;
	229	}
	230
	231	Bool_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
	232	{
	233	return GetDistanceFromVertexXY( vtx, 0, val, err );
	234	}
	235
	236
	237	Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx, Double_t &val, Double_t &err ) const
e7b09c95	238	{
	239	//* Calculate distance from vertex [cm] in XY-plane
	240
446ce366	241	return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
	242	}
	243
	244	Bool_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx, Double_t &val, Double_t &err ) const
	245	{
	246	//* Calculate distance from vertex [cm] in XY-plane
	247
	248	return GetDistanceFromVertexXY( AliKFParticle(Vtx), val, err );
	249	}
	250
	251	Double_t AliKFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
	252	{
	253	//* Calculate distance from vertex [cm] in XY-plane
	254	Double_t val, err;
	255	GetDistanceFromVertexXY( vtx, 0, val, err );
	256	return val;
e7b09c95	257	}
	258
	259	Double_t AliKFParticle::GetDistanceFromVertexXY( const AliKFParticle &Vtx ) const
	260	{
	261	//* Calculate distance from vertex [cm] in XY-plane
	262
	263	return GetDistanceFromVertexXY( Vtx.fP );
	264	}
	265
706952f5	266	Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const
e7b09c95	267	{
	268	//* Calculate distance from vertex [cm] in XY-plane
	269
	270	return GetDistanceFromVertexXY( AliKFParticle(Vtx).fP );
	271	}
	272
	273	Double_t AliKFParticle::GetDistanceFromParticleXY( const AliKFParticle &p ) const
	274	{
	275	//* Calculate distance to other particle [cm]
	276
	277	Double_t dS, dS1;
	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(dxdx+dydy);
	285	}
	286
	287	Double_t AliKFParticle::GetDeviationFromParticleXY( const AliKFParticle &p ) const
	288	{
	289	//* Calculate sqrt(Chi2/ndf) deviation from other particle
	290
	291	Double_t dS, dS1;
	292	GetDStoParticleXY( p, dS, dS1 );
	293	Double_t mP1[8], mC1[36];
	294	p.Transport( dS1, mP1, mC1 );
616ffc76	295
e7b09c95	296	Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
616ffc76	297
e7b09c95	298	Double_t sigmaS = .1+10.TMath::Sqrt( (d[0]d[0]+d[1]*d[1] )/
e7b09c95	299	(mP1[3]mP1[3]+mP1[4]mP1[4] ) );
616ffc76	300
e7b09c95	301	Double_t h[2] = { mP1[3]sigmaS, mP1[4]sigmaS };
616ffc76	302
e7b09c95	303	mC1[0] +=h[0]*h[0];
	304	mC1[1] +=h[1]*h[0];
	305	mC1[2] +=h[1]*h[1];
	306
	307	return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
	308	}
	309
	310
446ce366	311	Double_t AliKFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const
e7b09c95	312	{
	313	//* Calculate sqrt(Chi2/ndf) deviation from vertex
	314	//* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
	315
446ce366	316	Double_t val, err;
	317	Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
	318	if( problem \|\| err<1.e-20 ) return 1.e4;
	319	else return val/err;
e7b09c95	320	}
	321
	322
706952f5	323	Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
e7b09c95	324	{
	325	//* Calculate sqrt(Chi2/ndf) deviation from vertex
	326	//* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
	327
	328	return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
	329	}
	330
706952f5	331	Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const
e7b09c95	332	{
	333	//* Calculate sqrt(Chi2/ndf) deviation from vertex
	334	//* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
	335
	336	AliKFParticle v(Vtx);
	337	return GetDeviationFromVertexXY( v.fP, v.fC );
	338	}
	339
e7b09c95	340	Double_t AliKFParticle::GetAngle ( const AliKFParticle &p ) const
	341	{
	342	//* Calculate the opening angle between two particles
	343
	344	Double_t dS, dS1;
	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] );
	351	n*=n1;
	352	Double_t a = 0;
	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();
	356	return a;
	357	}
	358
	359	Double_t AliKFParticle::GetAngleXY( const AliKFParticle &p ) const
	360	{
	361	//* Calculate the opening angle between two particles in XY plane
	362
	363	Double_t dS, dS1;
	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] );
	370	n*=n1;
	371	Double_t a = 0;
	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();
	375	return a;
	376	}
	377
	378	Double_t AliKFParticle::GetAngleRZ( const AliKFParticle &p ) const
	379	{
	380	//* Calculate the opening angle between two particles in RZ plane
616ffc76	381
e7b09c95	382	Double_t dS, dS1;
	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( nrnr + mP[5]mP[5] );
	390	Double_t n1= TMath::Sqrt( n1rn1r + mP1[5]mP1[5] );
	391	n*=n1;
	392	Double_t a = 0;
	393	if( n>1.e-8 ) a = ( nrn1r +mP[5]mP1[5])/n;
	394	if (TMath::Abs(a)<1.) a = TMath::ACos(a);
	395	else a = (a>=0) ?0 :TMath::Pi();
	396	return a;
616ffc76	397	}
706952f5	398
	399
	400	/*
	401
	402	#include "AliExternalTrackParam.h"
	403
	404	void AliKFParticle::GetDStoParticleALICE( const AliKFParticleBase &p,
	405	Double_t &DS, Double_t &DS1 )
	406	const
	407	{
	408	DS = DS1 = 0;
	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;
	413
	414	GetExternalTrackParam( *this, x1, a1, par1 );
	415	GetExternalTrackParam( p, x2, a2, par2 );
	416
	417	AliExternalTrackParam t1(x1,a1, par1, cov);
	418	AliExternalTrackParam t2(x2,a2, par2, cov);
	419
	420	Double_t xe1=0, xe2=0;
	421	t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
	422	t1.PropagateTo( xe1, -GetFieldAlice() );
	423	t2.PropagateTo( xe2, -GetFieldAlice() );
	424
	425	Double_t xyz1[3], xyz2[3];
	426	t1.GetXYZ( xyz1 );
	427	t2.GetXYZ( xyz2 );
	428
	429	DS = GetDStoPoint( xyz1 );
	430	DS1 = p.GetDStoPoint( xyz2 );
	431
	432	return;
	433	}
	434	*/
91b25235	435
	436	// * Pseudo Proper Time of decay = (rpt) / \|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();
	443
	444	if ( timeErr2 ) {
	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,
	448	// py*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),
	451	// -px*m/pt^2,
	452	// -py*m/pt^2 }
	453	const Double_t f0 = Px()*mipt2;
	454	const Double_t f1 = Py()*mipt2;
	455	const Double_t mipt2derivative = mipt2(1-2Px()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;
	460
5f653c6e	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);
91b25235	474
91b25235	475	*timeErr2 =
5f653c6e	476	f5mC55f5 +
	477	f5mC54f4 +
	478	f4mC44f4 +
	479	f3mC33f3 +
	480	f3mC32f2 +
	481	f3mC31f1 +
	482	f3mC30f0 +
	483	f2mC22f2 +
	484	f2mC21f1 +
	485	f2mC20f0 +
	486	f1mC11f1 +
	487	f1mC10f0 +
	488	f0mC00f0;
91b25235	489	}
	490	return ( dxPx() + dyPy() )*mipt2;
	491	}