StRoot: StarRoot/KFParticle.cxx Source File

00001 //----------------------------------------------------------------------------
00002 // Implementation of the KFParticle class
00003 // .
00004 // @author  S.Gorbunov, I.Kisel
00005 // @version 1.0
00006 // @since   13.05.07
00007 // 
00008 // Class to reconstruct and store the decayed particle parameters.
00009 // The method is described in CBM-SOFT note 2007-003, 
00010 // ``Reconstruction of decayed particles based on the Kalman filter'', 
00011 // http://www.gsi.de/documents/DOC-2007-May-14-1.pdf
00012 //
00013 // This class is ALICE interface to general mathematics in KFParticleCore
00014 // 
00015 //  -= Copyright &copy ALICE HLT Group =-
00016 //____________________________________________________________________________
00017 
00018 
00019 #include "KFParticle.h"
00020 #include "TDatabasePDG.h"
00021 #include "TParticlePDG.h"
00022 #include "MTrack.h"
00023 #include "MVertex.h"
00024 ClassImp(KFParticle);
00025 
00026 
00027 Double_t KFParticle::fgBz = -5.;  //* Bz compoment of the magnetic field
00028 
00029 void KFParticle::Create( const Double_t Param[], const Double_t Cov[], Int_t Charge, Int_t PID )
00030 {
00031   // Constructor from "cartesian" track, PID hypothesis should be provided
00032   //
00033   // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
00034   // Cov [21] = lower-triangular part of the covariance matrix:
00035   //
00036   //                (  0  .  .  .  .  . )
00037   //                (  1  2  .  .  .  . )
00038   //  Cov. matrix = (  3  4  5  .  .  . ) - numbering of covariance elements in Cov[]
00039   //                (  6  7  8  9  .  . )
00040   //                ( 10 11 12 13 14  . )
00041   //                ( 15 16 17 18 19 20 )
00042   Double_t C[21];
00043   for( int i=0; i<21; i++ ) C[i] = Cov[i];
00044   
00045   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(PID);
00046   Double_t mass = (particlePDG) ? particlePDG->Mass() :0.13957;
00047   
00048   KFParticleBase::Initialize( Param, C, Charge, mass );
00049 }
00050 
00051 KFParticle::KFParticle( const MTrack &track, Int_t PID )
00052 {
00053   // Constructor from ALICE track, PID hypothesis should be provided
00054 
00055   track.XvYvZv(fP);
00056   track.PxPyPz(fP+3);
00057   fQ = track.Charge();
00058   track.GetCovarianceXYZPxPyPz( fC );
00059   Create(fP,fC,fQ,PID);
00060 }
00061 
00062 KFParticle::KFParticle( const MVertex &vertex )
00063 {
00064   // Constructor from ALICE vertex
00065 
00066   vertex.GetXYZ( fP );
00067   vertex.GetCovarianceMatrix( fC );  
00068   fChi2 = vertex.GetChi2();
00069   fNDF = 2*vertex.GetNContributors() - 3;
00070   fQ = 0;
00071   fAtProductionVertex = 0;
00072   fIsLinearized = 0;
00073   fSFromDecay = 0;
00074 }
00075 
00076 void KFParticle::GetExternalTrackParam( const KFParticleBase &p, Double_t &X, Double_t &Alpha, Double_t P[5] ) 
00077 {
00078   // Conversion to AliExternalTrackParam parameterization
00079 
00080   Double_t cosA = p.GetPx(), sinA = p.GetPy(); 
00081   Double_t pt = TMath::Sqrt(cosA*cosA + sinA*sinA);
00082   Double_t pti = 0;
00083   if( pt<1.e-4 ){
00084     cosA = 1;
00085     sinA = 0;
00086   } else {
00087     pti = 1./pt;
00088     cosA*=pti;
00089     sinA*=pti;
00090   }
00091   Alpha = TMath::ATan2(sinA,cosA);  
00092   X   = p.GetX()*cosA + p.GetY()*sinA;   
00093   P[0]= p.GetY()*cosA - p.GetX()*sinA;
00094   P[1]= p.GetZ();
00095   P[2]= 0;
00096   P[3]= p.GetPz()*pti;
00097   P[4]= p.GetQ()*pti;
00098 }
00099 
00100 Bool_t KFParticle::GetDistanceFromVertexXY( const Double_t vtx[], const Double_t Cv[], Double_t &val, Double_t &err ) const
00101 {
00102   //* Calculate DCA distance from vertex (transverse impact parameter) in XY
00103   //* v = [xy], Cv=[Cxx,Cxy,Cyy ]-covariance matrix
00104 
00105   Bool_t ret = 0;
00106   
00107   Double_t mP[8];
00108   Double_t mC[36];
00109   
00110   Transport( GetDStoPoint(vtx), mP, mC );  
00111 
00112   Double_t dx = mP[0] - vtx[0];
00113   Double_t dy = mP[1] - vtx[1];
00114   Double_t px = mP[3];
00115   Double_t py = mP[4];
00116   Double_t pt = TMath::Sqrt(px*px + py*py);
00117   Double_t ex=0, ey=0;
00118   if( pt<1.e-4 ){
00119     ret = 1;
00120     pt = 1.;
00121     val = 1.e4;
00122   } else{
00123     ex = px/pt;
00124     ey = py/pt;
00125     val = dy*ex - dx*ey;
00126   }
00127 
00128   Double_t h0 = -ey;
00129   Double_t h1 = ex;
00130   Double_t h3 = (dy*ey + dx*ex)*ey/pt;
00131   Double_t h4 = -(dy*ey + dx*ex)*ex/pt;
00132   
00133   err = 
00134     h0*(h0*GetCovariance(0,0) + h1*GetCovariance(0,1) + h3*GetCovariance(0,3) + h4*GetCovariance(0,4) ) +
00135     h1*(h0*GetCovariance(1,0) + h1*GetCovariance(1,1) + h3*GetCovariance(1,3) + h4*GetCovariance(1,4) ) +
00136     h3*(h0*GetCovariance(3,0) + h1*GetCovariance(3,1) + h3*GetCovariance(3,3) + h4*GetCovariance(3,4) ) +
00137     h4*(h0*GetCovariance(4,0) + h1*GetCovariance(4,1) + h3*GetCovariance(4,3) + h4*GetCovariance(4,4) );
00138 
00139   if( Cv ){
00140     err+= h0*(h0*Cv[0] + h1*Cv[1] ) + h1*(h0*Cv[1] + h1*Cv[2] ); 
00141   }
00142 
00143   err = TMath::Sqrt(TMath::Abs(err));
00144 
00145   return ret;
00146 }
00147 
00148 Bool_t KFParticle::GetDistanceFromVertexXY( const Double_t vtx[], Double_t &val, Double_t &err ) const
00149 {
00150   return GetDistanceFromVertexXY( vtx, 0, val, err );
00151 }
00152 
00153 
00154 Bool_t KFParticle::GetDistanceFromVertexXY( const KFParticle &Vtx, Double_t &val, Double_t &err ) const 
00155 {
00156   //* Calculate distance from vertex [cm] in XY-plane
00157 
00158   return GetDistanceFromVertexXY( Vtx.fP, Vtx.fC, val, err );
00159 }
00160 
00161 Bool_t KFParticle::GetDistanceFromVertexXY( const MVertex &Vtx, Double_t &val, Double_t &err ) const 
00162 {
00163   //* Calculate distance from vertex [cm] in XY-plane
00164 
00165   return GetDistanceFromVertexXY( KFParticle(Vtx), val, err );
00166 }
00167 
00168 Double_t KFParticle::GetDistanceFromVertexXY( const Double_t vtx[] ) const
00169 {
00170   //* Calculate distance from vertex [cm] in XY-plane
00171   Double_t val, err;
00172   GetDistanceFromVertexXY( vtx, 0, val, err );
00173   return val;
00174 }
00175 
00176 Double_t KFParticle::GetDistanceFromVertexXY( const KFParticle &Vtx ) const 
00177 {
00178   //* Calculate distance from vertex [cm] in XY-plane
00179 
00180   return GetDistanceFromVertexXY( Vtx.fP );
00181 }
00182 
00183 Double_t KFParticle::GetDistanceFromVertexXY( const MVertex &Vtx ) const 
00184 {
00185   //* Calculate distance from vertex [cm] in XY-plane
00186 
00187   return GetDistanceFromVertexXY( KFParticle(Vtx).fP );
00188 }
00189 
00190 Double_t KFParticle::GetDistanceFromParticleXY( const KFParticle &p ) const 
00191 {
00192   //* Calculate distance to other particle [cm]
00193 
00194   Double_t dS, dS1;
00195   GetDStoParticleXY( p, dS, dS1 );   
00196   Double_t mP[8], mC[36], mP1[8], mC1[36];
00197   Transport( dS, mP, mC ); 
00198   p.Transport( dS1, mP1, mC1 ); 
00199   Double_t dx = mP[0]-mP1[0]; 
00200   Double_t dy = mP[1]-mP1[1]; 
00201   return TMath::Sqrt(dx*dx+dy*dy);
00202 }
00203 
00204 Double_t KFParticle::GetDeviationFromParticleXY( const KFParticle &p ) const 
00205 {
00206   //* Calculate sqrt(Chi2/ndf) deviation from other particle
00207 
00208   Double_t dS, dS1;
00209   GetDStoParticleXY( p, dS, dS1 );   
00210   Double_t mP1[8], mC1[36];
00211   p.Transport( dS1, mP1, mC1 ); 
00212 
00213   Double_t d[2]={ fP[0]-mP1[0], fP[1]-mP1[1] };
00214 
00215   Double_t sigmaS = .1+10.*TMath::Sqrt( (d[0]*d[0]+d[1]*d[1] )/
00216                                         (mP1[3]*mP1[3]+mP1[4]*mP1[4] )  );
00217 
00218   Double_t h[2] = { mP1[3]*sigmaS, mP1[4]*sigmaS };       
00219   
00220   mC1[0] +=h[0]*h[0];
00221   mC1[1] +=h[1]*h[0]; 
00222   mC1[2] +=h[1]*h[1]; 
00223 
00224   return GetDeviationFromVertexXY( mP1, mC1 )*TMath::Sqrt(2./1.);
00225 }
00226 
00227 
00228 Double_t KFParticle::GetDeviationFromVertexXY( const Double_t vtx[], const Double_t Cv[] ) const 
00229 {
00230   //* Calculate sqrt(Chi2/ndf) deviation from vertex
00231   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
00232 
00233   Double_t val, err;
00234   Bool_t problem = GetDistanceFromVertexXY( vtx, Cv, val, err );
00235   if( problem || err<1.e-20 ) return 1.e4;
00236   else return val/err;
00237 }
00238 
00239 
00240 Double_t KFParticle::GetDeviationFromVertexXY( const KFParticle &Vtx ) const  
00241 {
00242   //* Calculate sqrt(Chi2/ndf) deviation from vertex
00243   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
00244 
00245   return GetDeviationFromVertexXY( Vtx.fP, Vtx.fC );
00246 }
00247 
00248 Double_t KFParticle::GetDeviationFromVertexXY( const MVertex &Vtx ) const 
00249 {
00250   //* Calculate sqrt(Chi2/ndf) deviation from vertex
00251   //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
00252 
00253   KFParticle v(Vtx);
00254   return GetDeviationFromVertexXY( v.fP, v.fC );
00255 }
00256 
00257 Double_t KFParticle::GetAngle  ( const KFParticle &p ) const 
00258 {
00259   //* Calculate the opening angle between two particles
00260 
00261   Double_t dS, dS1;
00262   GetDStoParticle( p, dS, dS1 );   
00263   Double_t mP[8], mC[36], mP1[8], mC1[36];
00264   Transport( dS, mP, mC ); 
00265   p.Transport( dS1, mP1, mC1 ); 
00266   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] + mP[5]*mP[5] );
00267   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] + mP1[5]*mP1[5] );
00268   n*=n1;
00269   Double_t a = 0;
00270   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] + mP[5]*mP1[5] )/n;
00271   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
00272   else a = (a>=0) ?0 :TMath::Pi();
00273   return a;
00274 }
00275 
00276 Double_t KFParticle::GetAngleXY( const KFParticle &p ) const 
00277 {
00278   //* Calculate the opening angle between two particles in XY plane
00279 
00280   Double_t dS, dS1;
00281   GetDStoParticleXY( p, dS, dS1 );   
00282   Double_t mP[8], mC[36], mP1[8], mC1[36];
00283   Transport( dS, mP, mC ); 
00284   p.Transport( dS1, mP1, mC1 ); 
00285   Double_t n = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
00286   Double_t n1= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4] );
00287   n*=n1;
00288   Double_t a = 0;
00289   if( n>1.e-8 ) a = ( mP[3]*mP1[3] + mP[4]*mP1[4] )/n;
00290   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
00291   else a = (a>=0) ?0 :TMath::Pi();
00292   return a;
00293 }
00294 
00295 Double_t KFParticle::GetAngleRZ( const KFParticle &p ) const 
00296 {
00297   //* Calculate the opening angle between two particles in RZ plane
00298 
00299   Double_t dS, dS1;
00300   GetDStoParticle( p, dS, dS1 );   
00301   Double_t mP[8], mC[36], mP1[8], mC1[36];
00302   Transport( dS, mP, mC ); 
00303   p.Transport( dS1, mP1, mC1 ); 
00304   Double_t nr = TMath::Sqrt( mP[3]*mP[3] + mP[4]*mP[4] );
00305   Double_t n1r= TMath::Sqrt( mP1[3]*mP1[3] + mP1[4]*mP1[4]  );
00306   Double_t n = TMath::Sqrt( nr*nr + mP[5]*mP[5] );
00307   Double_t n1= TMath::Sqrt( n1r*n1r + mP1[5]*mP1[5] );
00308   n*=n1;
00309   Double_t a = 0;
00310   if( n>1.e-8 ) a = ( nr*n1r +mP[5]*mP1[5])/n; 
00311   if (TMath::Abs(a)<1.) a = TMath::ACos(a);
00312   else a = (a>=0) ?0 :TMath::Pi();
00313   return a;
00314 }
00315 
00316 
00317 /*
00318 
00319 #include "AliExternalTrackParam.h"
00320 
00321 void KFParticle::GetDStoParticleALICE( const KFParticleBase &p, 
00322                                            Double_t &DS, Double_t &DS1 ) 
00323   const
00324 { 
00325   DS = DS1 = 0;   
00326   Double_t x1, a1, x2, a2;
00327   Double_t par1[5], par2[5], cov[15];
00328   for(int i=0; i<15; i++) cov[i] = 0;
00329   cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = .001;
00330 
00331   GetExternalTrackParam( *this, x1, a1, par1 );
00332   GetExternalTrackParam( p, x2, a2, par2 );
00333 
00334   AliExternalTrackParam t1(x1,a1, par1, cov);
00335   AliExternalTrackParam t2(x2,a2, par2, cov);
00336 
00337   Double_t xe1=0, xe2=0;
00338   t1.GetDCA( &t2, -GetFieldAlice(), xe1, xe2 );
00339   t1.PropagateTo( xe1, -GetFieldAlice() );
00340   t2.PropagateTo( xe2, -GetFieldAlice() );
00341 
00342   Double_t xyz1[3], xyz2[3];
00343   t1.GetXYZ( xyz1 );
00344   t2.GetXYZ( xyz2 );
00345   
00346   DS = GetDStoPoint( xyz1 );
00347   DS1 = p.GetDStoPoint( xyz2 );
00348 
00349   return;
00350 }
00351 */