]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliExternalTrackParam.cxx
Typo fixed.
[u/mrichter/AliRoot.git] / STEER / AliExternalTrackParam.cxx
index 2dc222ec7b0104cf32908fb0309572df6f14389b..e60bd30ef4cb64ce88c97c4db8682bac7b2dd264 100644 (file)
 // are implemented.
 // Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                            //
 ///////////////////////////////////////////////////////////////////////////////
+#include <TMatrixDSym.h>
 #include "AliExternalTrackParam.h"
 #include "AliESDVertex.h"
+#include "TPolyMarker3D.h"
 #include "AliLog.h"
 
 ClassImp(AliExternalTrackParam)
 
+Double32_t AliExternalTrackParam::fgMostProbablePt=kMostProbablePt;
 //_____________________________________________________________________________
 AliExternalTrackParam::AliExternalTrackParam() :
-  TObject(),
+  AliVParticle(),
   fX(0),
   fAlpha(0)
 {
@@ -46,7 +50,7 @@ AliExternalTrackParam::AliExternalTrackParam() :
 
 //_____________________________________________________________________________
 AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track):
-  TObject(track),
+  AliVParticle(track),
   fX(track.fX),
   fAlpha(track.fAlpha)
 {
@@ -57,11 +61,30 @@ AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track)
   for (Int_t i = 0; i < 15; i++) fC[i] = track.fC[i];
 }
 
+//_____________________________________________________________________________
+AliExternalTrackParam& AliExternalTrackParam::operator=(const AliExternalTrackParam &trkPar)
+{
+  //
+  // assignment operator
+  //
+  
+  if (this!=&trkPar) {
+    AliVParticle::operator=(trkPar);
+    fX = trkPar.fX;
+    fAlpha = trkPar.fAlpha;
+
+    for (Int_t i = 0; i < 5; i++) fP[i] = trkPar.fP[i];
+    for (Int_t i = 0; i < 15; i++) fC[i] = trkPar.fC[i];
+  }
+
+  return *this;
+}
+
 //_____________________________________________________________________________
 AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha, 
                                             const Double_t param[5], 
                                             const Double_t covar[15]) :
-  TObject(),
+  AliVParticle(),
   fX(x),
   fAlpha(alpha)
 {
@@ -94,6 +117,20 @@ void AliExternalTrackParam::Reset() {
   for (Int_t i = 0; i < 15; i++) fC[i] = 0;
 }
 
+//_____________________________________________________________________________
+void AliExternalTrackParam::AddCovariance(const Double_t c[15]) {
+  //
+  // Add "something" to the track covarince matrix.
+  // May be needed to account for unknown mis-calibration/mis-alignment
+  //
+    fC[0] +=c[0];
+    fC[1] +=c[1];  fC[2] +=c[2];
+    fC[3] +=c[3];  fC[4] +=c[4];  fC[5] +=c[5];
+    fC[6] +=c[6];  fC[7] +=c[7];  fC[8] +=c[8];  fC[9] +=c[9];
+    fC[10]+=c[10]; fC[11]+=c[11]; fC[12]+=c[12]; fC[13]+=c[13]; fC[14]+=c[14];
+}
+
+
 Double_t AliExternalTrackParam::GetP() const {
   //---------------------------------------------------------------------
   // This function returns the track momentum
@@ -178,9 +215,72 @@ Double_t AliExternalTrackParam::GetLinearD(Double_t xv,Double_t yv) const {
   return -d;
 }
 
-Bool_t AliExternalTrackParam::
-CorrectForMaterial(Double_t d,  Double_t x0, Double_t mass) {
+Bool_t AliExternalTrackParam::CorrectForMeanMaterial
+(Double_t xOverX0,  Double_t xTimesRho, Double_t mass, Bool_t anglecorr, 
+ Double_t (*Bethe)(Double_t)) {
+  //------------------------------------------------------------------
+  // This function corrects the track parameters for the crossed material.
+  // "xOverX0"   - X/X0, the thickness in units of the radiation length.
+  // "xTimesRho" - is the product length*density (g/cm^2). 
+  // "mass" - the mass of this particle (GeV/c^2).
+  //------------------------------------------------------------------
+  Double_t &fP2=fP[2];
+  Double_t &fP3=fP[3];
+  Double_t &fP4=fP[4];
+
+  Double_t &fC22=fC[5];
+  Double_t &fC33=fC[9];
+  Double_t &fC43=fC[13];
+  Double_t &fC44=fC[14];
+
+  //Apply angle correction, if requested
+  if(anglecorr) {
+    Double_t angle=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
+    xOverX0 *=angle;
+    xTimesRho *=angle;
+  } 
+
+  Double_t p=GetP();
+  Double_t p2=p*p;
+  Double_t beta2=p2/(p2 + mass*mass);
+
+  //Multiple scattering******************
+  if (xOverX0 != 0) {
+     Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(xOverX0);
+     if(theta2>TMath::Pi()*TMath::Pi()) return kFALSE;
+     //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
+     fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
+     fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
+     fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
+     fC44 += theta2*fP3*fP4*fP3*fP4;
+  }
+
+  //Energy losses************************
+  if ((xTimesRho != 0.) && (beta2 < 1.)) {
+     Double_t dE=Bethe(beta2)*xTimesRho;
+     Double_t e=TMath::Sqrt(p2 + mass*mass);
+     if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
+     fP4*=(1.- e/p2*dE);
+     if (TMath::Abs(fP4)>100.) return kFALSE; // Do not track below 10 MeV/c
+
+
+     // Approximate energy loss fluctuation (M.Ivanov)
+     const Double_t knst=0.07; // To be tuned.  
+     Double_t sigmadE=knst*TMath::Sqrt(TMath::Abs(dE)); 
+     fC44+=((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4)); 
+  }
+
+  return kTRUE;
+}
+
+
+Bool_t AliExternalTrackParam::CorrectForMaterial
+(Double_t d,  Double_t x0, Double_t mass, Double_t (*Bethe)(Double_t)) {
   //------------------------------------------------------------------
+  //                    Deprecated function !   
+  //       Better use CorrectForMeanMaterial instead of it.
+  //
   // This function corrects the track parameters for the crossed material
   // "d"    - the thickness (fraction of the radiation length)
   // "x0"   - the radiation length (g/cm^2) 
@@ -203,6 +303,7 @@ CorrectForMaterial(Double_t d,  Double_t x0, Double_t mass) {
   //Multiple scattering******************
   if (d!=0) {
      Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
+     if(theta2>TMath::Pi()*TMath::Pi()) return kFALSE;
      //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
      fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
      fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
@@ -213,16 +314,35 @@ CorrectForMaterial(Double_t d,  Double_t x0, Double_t mass) {
   //Energy losses************************
   if (x0!=0. && beta2<1) {
      d*=x0;
-     Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
-     if (beta2/(1-beta2)>3.5*3.5)
-       dE=0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2)*d;
-
-     fP4*=(1.- TMath::Sqrt(p2 + mass*mass)/p2*dE);
+     Double_t dE=Bethe(beta2)*d;
+     Double_t e=TMath::Sqrt(p2 + mass*mass);
+     if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
+     fP4*=(1.- e/p2*dE);
+
+     // Approximate energy loss fluctuation (M.Ivanov)
+     const Double_t knst=0.07; // To be tuned.  
+     Double_t sigmadE=knst*TMath::Sqrt(TMath::Abs(dE)); 
+     fC44+=((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4)); 
   }
 
   return kTRUE;
 }
 
+Double_t ApproximateBetheBloch(Double_t beta2) {
+  //------------------------------------------------------------------
+  // This is an approximation of the Bethe-Bloch formula with 
+  // the density effect taken into account at beta*gamma > 3.5
+  // (the approximation is reasonable only for solid materials) 
+  //------------------------------------------------------------------
+  if (beta2 >= 1) return kVeryBig;
+
+  if (beta2/(1-beta2)>3.5*3.5)
+     return 0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2);
+
+  return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2);
+}
+
 Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
   //------------------------------------------------------------------
   // Transform this track to the local coord. system rotated
@@ -364,7 +484,7 @@ Double_t p[3], Double_t bz) const {
   //+++++++++++++++++++++++++++++++++++++++++    
   GetXYZ(x);
     
-  if (TMath::Abs(Get1Pt()) < kAlmost0){ //straight-line tracks
+  if (OneOverPt() < kAlmost0 || TMath::Abs(bz) < kAlmost0Field ){ //straight-line tracks
      Double_t unit[3]; GetDirection(unit);
      x[0]+=unit[0]*len;   
      x[1]+=unit[1]*len;   
@@ -440,6 +560,148 @@ AliExternalTrackParam::GetPredictedChi2(Double_t p[2],Double_t cov[3]) const {
   return (d*szz*d - 2*d*sdz*z + z*sdd*z)/det;
 }
 
+Double_t AliExternalTrackParam::
+GetPredictedChi2(Double_t p[3],Double_t covyz[3],Double_t covxyz[3]) const {
+  //----------------------------------------------------------------
+  // Estimate the chi2 of the 3D space point "p" and
+  // the full covariance matrix "covyz" and "covxyz"
+  //
+  // Cov(x,x) ... :   covxyz[0]
+  // Cov(y,x) ... :   covxyz[1]  covyz[0]
+  // Cov(z,x) ... :   covxyz[2]  covyz[1]  covyz[2]
+  //----------------------------------------------------------------
+
+  Double_t res[3] = {
+    GetX() - p[0],
+    GetY() - p[1],
+    GetZ() - p[2]
+  };
+
+  Double_t f=GetSnp();
+  if (TMath::Abs(f) >= kAlmost1) return kVeryBig;
+  Double_t r=TMath::Sqrt(1.- f*f);
+  Double_t a=f/r, b=GetTgl()/r;
+
+  Double_t s2=333.*333.;  //something reasonably big (cm^2)
+  TMatrixDSym v(3);
+  v(0,0)=  s2;  v(0,1)=  a*s2;                 v(0,2)=  b*s2;;
+  v(1,0)=a*s2;  v(1,1)=a*a*s2 + GetSigmaY2();  v(1,2)=a*b*s2 + GetSigmaZY();
+  v(2,0)=b*s2;  v(2,1)=a*b*s2 + GetSigmaZY();  v(2,2)=b*b*s2 + GetSigmaZ2();
+
+  v(0,0)+=covxyz[0]; v(0,1)+=covxyz[1]; v(0,2)+=covxyz[2];
+  v(1,0)+=covxyz[1]; v(1,1)+=covyz[0];  v(1,2)+=covyz[1];
+  v(2,0)+=covxyz[2]; v(2,1)+=covyz[1];  v(2,2)+=covyz[2];
+
+  v.Invert();
+  if (!v.IsValid()) return kVeryBig;
+
+  Double_t chi2=0.;
+  for (Int_t i = 0; i < 3; i++)
+    for (Int_t j = 0; j < 3; j++) chi2 += res[i]*res[j]*v(i,j);
+
+  return chi2;  
+
+
+}
+
+Bool_t AliExternalTrackParam::
+PropagateTo(Double_t p[3],Double_t covyz[3],Double_t covxyz[3],Double_t bz) {
+  //----------------------------------------------------------------
+  // Propagate this track to the plane 
+  // the 3D space point "p" (with the covariance matrix "covyz" and "covxyz")
+  // belongs to.
+  // The magnetic field is "bz" (kG)
+  //
+  // The track curvature and the change of the covariance matrix
+  // of the track parameters are negleted !
+  // (So the "step" should be small compared with 1/curvature)
+  //----------------------------------------------------------------
+
+  Double_t f=GetSnp();
+  if (TMath::Abs(f) >= kAlmost1) return kFALSE;
+  Double_t r=TMath::Sqrt(1.- f*f);
+  Double_t a=f/r, b=GetTgl()/r;
+
+  Double_t s2=333.*333.;  //something reasonably big (cm^2)
+  TMatrixDSym tV(3);
+  tV(0,0)=  s2;  tV(0,1)=  a*s2;  tV(0,2)=  b*s2;
+  tV(1,0)=a*s2;  tV(1,1)=a*a*s2;  tV(1,2)=a*b*s2;
+  tV(2,0)=b*s2;  tV(2,1)=a*b*s2;  tV(2,2)=b*b*s2;
+
+  TMatrixDSym pV(3);
+  pV(0,0)=covxyz[0]; pV(0,1)=covxyz[1]; pV(0,2)=covxyz[2];
+  pV(1,0)=covxyz[1]; pV(1,1)=covyz[0];  pV(1,2)=covyz[1];
+  pV(2,0)=covxyz[2]; pV(2,1)=covyz[1];  pV(2,2)=covyz[2];
+
+  TMatrixDSym tpV(tV);
+  tpV+=pV;
+  tpV.Invert();
+  if (!tpV.IsValid()) return kFALSE;
+
+  TMatrixDSym pW(3),tW(3);
+  for (Int_t i=0; i<3; i++)
+    for (Int_t j=0; j<3; j++) {
+      pW(i,j)=tW(i,j)=0.;
+      for (Int_t k=0; k<3; k++) {
+       pW(i,j) += tV(i,k)*tpV(k,j);
+       tW(i,j) += pV(i,k)*tpV(k,j);
+      }
+    }
+
+  Double_t t[3] = {GetX(), GetY(), GetZ()};
+
+  Double_t x=0.;
+  for (Int_t i=0; i<3; i++) x += (tW(0,i)*t[i] + pW(0,i)*p[i]);  
+  Double_t crv=GetC(bz);
+  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
+  f += crv*(x-fX);
+  if (TMath::Abs(f) >= kAlmost1) return kFALSE;
+  fX=x;  
+
+  fP[0]=0.;
+  for (Int_t i=0; i<3; i++) fP[0] += (tW(1,i)*t[i] + pW(1,i)*p[i]);  
+  fP[1]=0.;
+  for (Int_t i=0; i<3; i++) fP[1] += (tW(2,i)*t[i] + pW(2,i)*p[i]);  
+
+  return kTRUE;  
+}
+
+Double_t *AliExternalTrackParam::GetResiduals(
+Double_t *p,Double_t *cov,Bool_t updated) const {
+  //------------------------------------------------------------------
+  // Returns the track residuals with the space point "p" having
+  // the covariance matrix "cov".
+  // If "updated" is kTRUE, the track parameters expected to be updated,
+  // otherwise they must be predicted.  
+  //------------------------------------------------------------------
+  static Double_t res[2];
+
+  Double_t r00=cov[0], r01=cov[1], r11=cov[2];
+  if (updated) {
+     r00-=fC[0]; r01-=fC[1]; r11-=fC[2];
+  } else {
+     r00+=fC[0]; r01+=fC[1]; r11+=fC[2];
+  }
+  Double_t det=r00*r11 - r01*r01;
+
+  if (TMath::Abs(det) < kAlmost0) return 0;
+
+  Double_t tmp=r00; r00=r11/det; r11=tmp/det;
+
+  if (r00 < 0.) return 0;
+  if (r11 < 0.) return 0;
+
+  Double_t dy = fP[0] - p[0];
+  Double_t dz = fP[1] - p[1];
+
+  res[0]=dy*TMath::Sqrt(r00);
+  res[1]=dz*TMath::Sqrt(r11);
+
+  return res;
+}
+
 Bool_t AliExternalTrackParam::Update(Double_t p[2], Double_t cov[3]) {
   //------------------------------------------------------------------
   // Update the track parameters with the space point "p" having
@@ -598,10 +860,10 @@ Double_t b, Double_t &xthis, Double_t &xp) const {
      if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
      if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
         if ((gt1*gt1+gt2*gt2) > 1.e-4/dy2/dy2) 
-         AliWarning(" stopped at not a stationary point !");
+         AliDebug(1," stopped at not a stationary point !");
         Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det);
         if (lmb < 0.) 
-         AliWarning(" stopped at not a minimum !");
+         AliDebug(1," stopped at not a minimum !");
         break;
      }
 
@@ -614,7 +876,7 @@ Double_t b, Double_t &xthis, Double_t &xp) const {
        if (dd<dm) break;
         dt1*=0.5; dt2*=0.5;
         if (div>512) {
-           AliWarning(" overshoot !"); break;
+         AliDebug(1," overshoot !"); break;
         }   
      }
      dm=dd;
@@ -624,7 +886,7 @@ Double_t b, Double_t &xthis, Double_t &xp) const {
 
   }
 
-  if (max<=0) AliWarning(" too many iterations !");
+  if (max<=0) AliDebug(1," too many iterations !");
 
   Double_t cs=TMath::Cos(GetAlpha());
   Double_t sn=TMath::Sin(GetAlpha());
@@ -661,11 +923,10 @@ PropagateToDCA(AliExternalTrackParam *p, Double_t b) {
 }
 
 
-
-
-Bool_t AliExternalTrackParam::PropagateToDCA(const AliESDVertex *vtx, Double_t b, Double_t maxd){
+Bool_t AliExternalTrackParam::PropagateToDCA(const AliESDVertex *vtx, 
+Double_t b, Double_t maxd, Double_t dz[2], Double_t covar[3]) {
   //
-  // Try to relate this track to the vertex "vtx", 
+  // Propagate this track to the DCA to vertex "vtx", 
   // if the (rough) transverse impact parameter is not bigger then "maxd". 
   //            Magnetic field is "b" (kG).
   //
@@ -679,7 +940,7 @@ Bool_t AliExternalTrackParam::PropagateToDCA(const AliESDVertex *vtx, Double_t b
   Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
   Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
   Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
-  Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs;
+  Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv();
   x-=xv; y-=yv;
 
   //Estimate the impact parameter neglecting the track curvature
@@ -687,65 +948,37 @@ Bool_t AliExternalTrackParam::PropagateToDCA(const AliESDVertex *vtx, Double_t b
   if (d > maxd) return kFALSE; 
 
   //Propagate to the DCA
-  Double_t crv=0.299792458e-3*b*GetParameter()[4];
+  Double_t crv=kB2C*b*GetParameter()[4];
+  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
+
   Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
   sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv); cs=TMath::Sqrt(1.- sn*sn);
+  if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
+  else cs=1.;
 
   x = xv*cs + yv*sn;
   yv=-xv*sn + yv*cs; xv=x;
 
   if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
-  return kTRUE;
-}
-
 
+  if (dz==0) return kTRUE;
+  dz[0] = GetParameter()[0] - yv;
+  dz[1] = GetParameter()[1] - zv;
+  
+  if (covar==0) return kTRUE;
+  Double_t cov[6]; vtx->GetCovMatrix(cov);
 
-
-Bool_t Local2GlobalMomentum(Double_t p[3],Double_t alpha) {
-  //----------------------------------------------------------------
-  // This function performs local->global transformation of the
-  // track momentum.
-  // When called, the arguments are:
-  //    p[0] = 1/pt of the track;
-  //    p[1] = sine of local azim. angle of the track momentum;
-  //    p[2] = tangent of the track momentum dip angle;
-  //   alpha - rotation angle. 
-  // The result is returned as:
-  //    p[0] = px
-  //    p[1] = py
-  //    p[2] = pz
-  // Results for (nearly) straight tracks are meaningless !
-  //----------------------------------------------------------------
-  if (TMath::Abs(p[0])<=kAlmost0) return kFALSE;
-  if (TMath::Abs(p[1])> kAlmost1) return kFALSE;
-
-  Double_t pt=1./TMath::Abs(p[0]);
-  Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
-  Double_t r=TMath::Sqrt(1 - p[1]*p[1]);
-  p[0]=pt*(r*cs - p[1]*sn); p[1]=pt*(p[1]*cs + r*sn); p[2]=pt*p[2];
+  //***** Improvements by A.Dainese
+  alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
+  Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
+  covar[0] = GetCovariance()[0] + s2ylocvtx;   // neglecting correlations
+  covar[1] = GetCovariance()[1];               // between (x,y) and z
+  covar[2] = GetCovariance()[2] + cov[5];      // in vertex's covariance matrix
+  //*****
 
   return kTRUE;
 }
 
-Bool_t Local2GlobalPosition(Double_t r[3],Double_t alpha) {
-  //----------------------------------------------------------------
-  // This function performs local->global transformation of the
-  // track position.
-  // When called, the arguments are:
-  //    r[0] = local x
-  //    r[1] = local y
-  //    r[2] = local z
-  //   alpha - rotation angle. 
-  // The result is returned as:
-  //    r[0] = global x
-  //    r[1] = global y
-  //    r[2] = global z
-  //----------------------------------------------------------------
-  Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha), x=r[0];
-  r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
-
-  return kTRUE;
-}
 
 void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
   //----------------------------------------------------------------
@@ -754,14 +987,14 @@ void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
   //----------------------------------------------------------------
   Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
   Double_t snp=fP[2];
-  Double_t csp =TMath::Sqrt(1.- snp*snp);
+  Double_t csp =TMath::Sqrt((1.- snp)*(1.+snp));
   Double_t norm=TMath::Sqrt(1.+ fP[3]*fP[3]);
   d[0]=(csp*cs - snp*sn)/norm; 
   d[1]=(snp*cs + csp*sn)/norm; 
   d[2]=fP[3]/norm;
 }
 
-Bool_t AliExternalTrackParam::GetPxPyPz(Double_t *p) const {
+Bool_t AliExternalTrackParam::GetPxPyPz(Double_t p[3]) const {
   //---------------------------------------------------------------------
   // This function returns the global track momentum components
   // Results for (nearly) straight tracks are meaningless !
@@ -770,6 +1003,127 @@ Bool_t AliExternalTrackParam::GetPxPyPz(Double_t *p) const {
   return Local2GlobalMomentum(p,fAlpha);
 }
 
+Double_t AliExternalTrackParam::Px() const {
+  //---------------------------------------------------------------------
+  // Returns x-component of momentum
+  // Result for (nearly) straight tracks is meaningless !
+  //---------------------------------------------------------------------
+
+  Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
+  GetPxPyPz(p);
+
+  return p[0];
+}
+
+Double_t AliExternalTrackParam::Py() const {
+  //---------------------------------------------------------------------
+  // Returns y-component of momentum
+  // Result for (nearly) straight tracks is meaningless !
+  //---------------------------------------------------------------------
+
+  Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
+  GetPxPyPz(p);
+
+  return p[1];
+}
+
+Double_t AliExternalTrackParam::Pz() const {
+  //---------------------------------------------------------------------
+  // Returns z-component of momentum
+  // Result for (nearly) straight tracks is meaningless !
+  //---------------------------------------------------------------------
+
+  Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
+  GetPxPyPz(p);
+
+  return p[2];
+}
+
+Double_t AliExternalTrackParam::Xv() const {
+  //---------------------------------------------------------------------
+  // Returns x-component of first track point
+  //---------------------------------------------------------------------
+
+  Double_t r[3]={0.,0.,0.};
+  GetXYZ(r);
+
+  return r[0];
+}
+
+Double_t AliExternalTrackParam::Yv() const {
+  //---------------------------------------------------------------------
+  // Returns y-component of first track point
+  //---------------------------------------------------------------------
+
+  Double_t r[3]={0.,0.,0.};
+  GetXYZ(r);
+
+  return r[1];
+}
+
+Double_t AliExternalTrackParam::Zv() const {
+  //---------------------------------------------------------------------
+  // Returns z-component of first track point
+  //---------------------------------------------------------------------
+
+  Double_t r[3]={0.,0.,0.};
+  GetXYZ(r);
+
+  return r[2];
+}
+
+Double_t AliExternalTrackParam::Theta() const {
+  // return theta angle of momentum
+
+  return 0.5*TMath::Pi() - TMath::ATan(fP[3]);
+}
+
+Double_t AliExternalTrackParam::Phi() const {
+  //---------------------------------------------------------------------
+  // Returns the azimuthal angle of momentum
+  // 0 <= phi < 2*pi
+  //---------------------------------------------------------------------
+
+  Double_t phi=TMath::ASin(fP[2]) + fAlpha;
+  if (phi<0.) phi+=2.*TMath::Pi();
+  else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi();
+  return phi;
+}
+
+Double_t AliExternalTrackParam::M() const {
+  // return particle mass
+
+  // No mass information available so far.
+  // Redifine in derived class!
+
+  return -999.;
+}
+
+Double_t AliExternalTrackParam::E() const {
+  // return particle energy
+
+  // No PID information available so far.
+  // Redifine in derived class!
+
+  return -999.;
+}
+
+Double_t AliExternalTrackParam::Eta() const { 
+  // return pseudorapidity
+
+  return -TMath::Log(TMath::Tan(0.5 * Theta())); 
+}
+
+Double_t AliExternalTrackParam::Y() const {
+  // return rapidity
+
+  // No PID information available so far.
+  // Redifine in derived class!
+
+  return -999.;
+}
+
 Bool_t AliExternalTrackParam::GetXYZ(Double_t *r) const {
   //---------------------------------------------------------------------
   // This function returns the global track position
@@ -801,13 +1155,17 @@ Bool_t AliExternalTrackParam::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
   }
   Double_t pt=1./TMath::Abs(fP[4]);
   Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
-  Double_t r=TMath::Sqrt(1-fP[2]*fP[2]);
+  Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
 
   Double_t m00=-sn, m10=cs;
   Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
   Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
   Double_t m35=pt, m45=-pt*pt*fP[3];
 
+  m43*=GetSign();
+  m44*=GetSign();
+  m45*=GetSign();
+
   cv[0 ] = fC[0]*m00*m00;
   cv[1 ] = fC[0]*m00*m10; 
   cv[2 ] = fC[0]*m10*m10;
@@ -934,3 +1292,77 @@ Double_t AliExternalTrackParam::GetSnpAt(Double_t x,Double_t b) const {
   Double_t res = fP[2]+dx*crv;
   return res;
 }
+
+Bool_t AliExternalTrackParam::GetDistance(AliExternalTrackParam *param2, Double_t x, Double_t dist[3], Double_t bz){
+  //------------------------------------------------------------------------
+  // Get the distance between two tracks at the local position x 
+  // working in the local frame of this track.
+  // Origin :   Marian.Ivanov@cern.ch
+  //-----------------------------------------------------------------------
+  Double_t xyz[3];
+  Double_t xyz2[3];
+  xyz[0]=x;
+  if (!GetYAt(x,bz,xyz[1])) return kFALSE;
+  if (!GetZAt(x,bz,xyz[2])) return kFALSE;
+  //  
+  //
+  if (TMath::Abs(GetAlpha()-param2->GetAlpha())<kAlmost0){
+    xyz2[0]=x;
+    if (!param2->GetYAt(x,bz,xyz2[1])) return kFALSE;
+    if (!param2->GetZAt(x,bz,xyz2[2])) return kFALSE;
+  }else{
+    //
+    Double_t xyz1[3];
+    Double_t dfi = param2->GetAlpha()-GetAlpha();
+    Double_t ca = TMath::Cos(dfi), sa = TMath::Sin(dfi);
+    xyz2[0] =  xyz[0]*ca+xyz[1]*sa;
+    xyz2[1] = -xyz[0]*sa+xyz[1]*ca;
+    //
+    xyz1[0]=xyz2[0];
+    if (!param2->GetYAt(xyz2[0],bz,xyz1[1])) return kFALSE;
+    if (!param2->GetZAt(xyz2[0],bz,xyz1[2])) return kFALSE;
+    //
+    xyz2[0] =  xyz1[0]*ca-xyz1[1]*sa;
+    xyz2[1] = +xyz1[0]*sa+xyz1[1]*ca;
+    xyz2[2] = xyz1[2];
+  }
+  dist[0] = xyz[0]-xyz2[0];
+  dist[1] = xyz[1]-xyz2[1];
+  dist[2] = xyz[2]-xyz2[2];
+
+  return kTRUE;
+}
+
+
+//
+// Draw functionality.
+// Origin: Marian Ivanov, Marian.Ivanov@cern.ch
+//
+
+void  AliExternalTrackParam::DrawTrack(Float_t magf, Float_t minR, Float_t maxR, Float_t stepR){
+  //
+  // Draw track line
+  //
+  if (minR>maxR) return ;
+  if (stepR<=0) return ;
+  Int_t npoints = TMath::Nint((maxR-minR)/stepR)+1;
+  if (npoints<1) return;
+  TPolyMarker3D *polymarker = new TPolyMarker3D(npoints);
+  FillPolymarker(polymarker, magf,minR,maxR,stepR);
+  polymarker->Draw();
+}
+
+//
+void AliExternalTrackParam::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
+  //
+  // Fill points in the polymarker
+  //
+  Int_t counter=0;
+  for (Double_t r=minR; r<maxR; r+=stepR){
+    Double_t point[3];
+    GetXYZAt(r,magF,point);
+    pol->SetPoint(counter,point[0],point[1], point[2]);
+    printf("xyz\t%f\t%f\t%f\n",point[0], point[1],point[2]);
+    counter++;
+  }
+}