[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 "AliExternalTrackParam.h"
//#include "TMath.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 )


  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
  Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
  
  AliKFParticleBase::Initialize( Param, Cov, Charge, mass );
}


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

  track.GetXYZ(fP);
  track.GetPxPyPz(fP+3);
  fQ = (track.GetSigned1Pt() >0 ) ?1 :-1;
  track.GetCovarianceXYZPxPyPz( fC );
  Create(fP,fC,fQ,PID);
}

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

  vertex.GetXYZ( fP );
  vertex.GetCovMatrix( 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 AliESDVertex &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 AliESDVertex &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;
}
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"
	22	#include "AliExternalTrackParam.h"
	23	//#include "TMath.h"
	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 )
	42
	43
f826d409	44	TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
f826d409	45	Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
e7b09c95	46
	47	AliKFParticleBase::Initialize( Param, Cov, Charge, mass );
	48	}
616ffc76	49
f826d409	50
e7b09c95	51	AliKFParticle::AliKFParticle( const AliExternalTrackParam &track, Int_t PID )
	52	{
	53	// Constructor from ALICE track, PID hypothesis should be provided
616ffc76	54
e7b09c95	55	track.GetXYZ(fP);
	56	track.GetPxPyPz(fP+3);
	57	fQ = (track.GetSigned1Pt() >0 ) ?1 :-1;
f826d409	58	track.GetCovarianceXYZPxPyPz( fC );
7972a8db	59	Create(fP,fC,fQ,PID);
f826d409	60	}
f826d409	61
4bbc290d	62	AliKFParticle::AliKFParticle( const AliESDVertex &vertex )
f826d409	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	}
616ffc76	75
	76
	77	void AliKFParticle::GetExternalTrackParam( const AliKFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] )
	78	{
e7b09c95	79	// Conversion to AliExternalTrackParam parameterization
e7b09c95	80
616ffc76	81	Double_t cosA = p.GetPx(), sinA = p.GetPy();
	82	Double_t pt = TMath::Sqrt(cosAcosA + sinAsinA);
	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
616ffc76	102
e7b09c95	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(dxdx+dydy);
	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 );
616ffc76	151
e7b09c95	152	Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
616ffc76	153
e7b09c95	154	Double_t sigmaS = .1+10.TMath::Sqrt( (d[0]d[0]+d[1]*d[1] )/
e7b09c95	155	(mP1[3]mP1[3]+mP1[4]mP1[4] ) );
616ffc76	156
e7b09c95	157	Double_t h[2] = { mP1[3]sigmaS, mP1[4]sigmaS };
616ffc76	158
e7b09c95	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
616ffc76	262
e7b09c95	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( nrnr + mP[5]mP[5] );
	271	Double_t n1= TMath::Sqrt( n1rn1r + mP1[5]mP1[5] );
	272	n*=n1;
	273	Double_t a = 0;
	274	if( n>1.e-8 ) a = ( nrn1r +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;
616ffc76	278	}