[u/mrichter/AliRoot.git] / STEER / 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"

ClassImp(AliKFParticle)

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

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 );
}


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

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

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;
}


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

  Double_t mP[8], mC[36];
  Transport( GetDStoPoint(vtx), mP, mC );
  Double_t d[2]={ vtx[0]-mP[0], vtx[1]-mP[1] };
  Double_t dist =  TMath::Sqrt( d[0]*d[0]+d[1]*d[1] );
  Double_t sign = d[0]*mP[3] - d[1]*mP[4];  
  return (sign>=0) ?dist :-dist;
}

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 v[], 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 mP[8];
  Double_t mC[36];
  
  Transport( GetDStoPoint(v), mP, mC );  

  Double_t d[2]={ v[0]-mP[0], v[1]-mP[1] };

  Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
					(mP[3]*mP[3]+mP[4]*mP[4] )  );
   
  Double_t h[2] = { mP[3]*sigmaS, mP[4]*sigmaS };       
  
  Double_t mSi[3] = { mC[0] +h[0]*h[0], 
		      mC[1] +h[1]*h[0], mC[2] +h[1]*h[1] };

  if( Cv ){
    mSi[0]+=Cv[0];
    mSi[1]+=Cv[1];
    mSi[2]+=Cv[2];
  }
  Double_t s = ( mSi[0]*mSi[2] - mSi[1]*mSi[1] );
  s = ( s > 1.E-20 )  ?1./s :0;	  
  
  Double_t mS[3] = { mSi[2], 
		     -mSi[1], mSi[0] };      

  return TMath::Sqrt( TMath::Abs(s*( ( mS[0]*d[0] + mS[1]*d[1] )*d[0]
				     +(mS[1]*d[0] + mS[2]*d[1] )*d[1] ))/1);
}


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;
}
*/
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
effa6338	25	ClassImp(AliKFParticle)
f826d409	26
616ffc76	27	Double_t AliKFParticle::fgBz = -5.; //* Bz compoment of the magnetic field
f826d409	28
7972a8db	29	void AliKFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
f826d409	30	{
e7b09c95	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 )
706952f5	42	Double_t C[21];
	43	for( int i=0; i<21; i++ ) C[i] = Cov[i];
	44
f826d409	45	TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
f826d409	46	Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
e7b09c95	47
706952f5	48	AliKFParticleBase::Initialize( Param, C, Charge, mass );
e7b09c95	49	}
616ffc76	50
f826d409	51
706952f5	52	AliKFParticle::AliKFParticle( const AliVTrack &track, Int_t PID )
e7b09c95	53	{
e7b09c95	54	// Constructor from ALICE track, PID hypothesis should be provided
616ffc76	55
706952f5	56	track.XvYvZv(fP);
	57	track.PxPyPz(fP+3);
	58	fQ = track.Charge();
f826d409	59	track.GetCovarianceXYZPxPyPz( fC );
7972a8db	60	Create(fP,fC,fQ,PID);
f826d409	61	}
f826d409	62
706952f5	63	AliKFParticle::AliKFParticle( const AliVVertex &vertex )
f826d409	64	{
	65	// Constructor from ALICE vertex
	66
	67	vertex.GetXYZ( fP );
706952f5	68	vertex.GetCovarianceMatrix( fC );
f826d409	69	fChi2 = vertex.GetChi2();
	70	fNDF = 2*vertex.GetNContributors() - 3;
	71	fQ = 0;
	72	fAtProductionVertex = 0;
	73	fIsLinearized = 0;
	74	fSFromDecay = 0;
	75	}
616ffc76	76
	77
	78	void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
	79	{
e7b09c95	80	// Conversion to AliExternalTrackParam parameterization
e7b09c95	81
616ffc76	82	Double_t cosA = p.GetPx(), sinA = p.GetPy();
	83	Double_t pt = TMath::Sqrt(cosAcosA + sinAsinA);
	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
616ffc76	103
e7b09c95	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
706952f5	123	Double_t AliKFParticle::GetDistanceFromVertexXY( const AliVVertex &Vtx ) const
e7b09c95	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(dxdx+dydy);
	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 );
616ffc76	152
e7b09c95	153	Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
616ffc76	154
e7b09c95	155	Double_t sigmaS = .1+10.TMath::Sqrt( (d[0]d[0]+d[1]*d[1] )/
e7b09c95	156	(mP1[3]mP1[3]+mP1[4]mP1[4] ) );
616ffc76	157
e7b09c95	158	Double_t h[2] = { mP1[3]sigmaS, mP1[4]sigmaS };
616ffc76	159
e7b09c95	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
706952f5	204	Double_t AliKFParticle::GetDeviationFromVertexXY( const AliKFParticle &Vtx ) const
e7b09c95	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
706952f5	212	Double_t AliKFParticle::GetDeviationFromVertexXY( const AliVVertex &Vtx ) const
e7b09c95	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
616ffc76	263
e7b09c95	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( nrnr + mP[5]mP[5] );
	272	Double_t n1= TMath::Sqrt( n1rn1r + mP1[5]mP1[5] );
	273	n*=n1;
	274	Double_t a = 0;
	275	if( n>1.e-8 ) a = ( nrn1r +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;
616ffc76	279	}
706952f5	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	*/