STEER/AliTrackerBase.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliTrackerBase.cxx 38069 2009-12-24 16:56:18Z belikov $ */
  17
  18 //-------------------------------------------------------------------------
  19 //               Implementation of the AliTrackerBase class
  20 //                that is the base for the AliTracker class
  21 //                     Origin: Marian.Ivanov@cern.ch
  22 //-------------------------------------------------------------------------
  23 #include <TClass.h>
  24 #include <TMath.h>
  25 #include <TGeoManager.h>
  26
  27 #include "AliLog.h"
  28 #include "AliTrackerBase.h"
  29 #include "AliExternalTrackParam.h"
  30 #include "AliTrackPointArray.h"
  31
  32 extern TGeoManager *gGeoManager;
  33
  34 ClassImp(AliTrackerBase)
  35
  36 AliTrackerBase::AliTrackerBase():
  37   TObject(),
  38   fX(0),
  39   fY(0),
  40   fZ(0),
  41   fSigmaX(0.005),
  42   fSigmaY(0.005),
  43   fSigmaZ(0.010)
  44 {
  45   //--------------------------------------------------------------------
  46   // The default constructor.
  47   //--------------------------------------------------------------------
  48   if (!TGeoGlobalMagField::Instance()->GetField())
  49     AliWarning("Field map is not set.");
  50 }
  51
  52 //__________________________________________________________________________
  53 AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr):
  54   TObject(atr),
  55   fX(atr.fX),
  56   fY(atr.fY),
  57   fZ(atr.fZ),
  58   fSigmaX(atr.fSigmaX),
  59   fSigmaY(atr.fSigmaY),
  60   fSigmaZ(atr.fSigmaZ)
  61 {
  62   //--------------------------------------------------------------------
  63   // The default constructor.
  64   //--------------------------------------------------------------------
  65   if (!TGeoGlobalMagField::Instance()->GetField())
  66     AliWarning("Field map is not set.");
  67 }
  68
  69 //__________________________________________________________________________
  70 Double_t AliTrackerBase::GetBz()
  71 {
  72   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  73   if (!fld) return 0.5*kAlmost0Field;
  74   Double_t bz = fld->SolenoidField();
  75   return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
  76 }
  77
  78 //__________________________________________________________________________
  79 Double_t AliTrackerBase::GetBz(const Double_t *r) {
  80   //------------------------------------------------------------------
  81   // Returns Bz (kG) at the point "r" .
  82   //------------------------------------------------------------------
  83   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  84   if (!fld) return  0.5*kAlmost0Field;
  85   Double_t bz = fld->GetBz(r);
  86   return  TMath::Sign(0.5*kAlmost0Field,bz) + bz;
  87 }
  88
  89 //__________________________________________________________________________
  90 void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) {
  91   //------------------------------------------------------------------
  92   // Returns Bx, By and Bz (kG) at the point "r" .
  93   //------------------------------------------------------------------
  94   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  95   if (!fld) {
  96      b[0] = b[1] = 0.;
  97      b[2] = 0.5*kAlmost0Field;
  98      return;
  99   }
 100
 101   if (fld->IsUniform()) {
 102      b[0] = b[1] = 0.;
 103      b[2] = fld->SolenoidField();
 104   }  else {
 105      fld->Field(r,b);
 106   }
 107   b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]);
 108   return;
 109 }
 110
 111 Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
 112 {
 113   //
 114   // Calculate mean material budget and material properties between
 115   //    the points "start" and "end".
 116   //
 117   // "mparam" - parameters used for the energy and multiple scattering
 118   //  corrections:
 119   //
 120   // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
 121   // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
 122   // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
 123   // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
 124   // mparam[4] - length: sum(x_i) [cm]
 125   // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
 126   // mparam[6] - number of boundary crosses
 127   //
 128   //  Origin:  Marian Ivanov, Marian.Ivanov@cern.ch
 129   //
 130   //  Corrections and improvements by
 131   //        Andrea Dainese, Andrea.Dainese@lnl.infn.it,
 132   //        Andrei Gheata,  Andrei.Gheata@cern.ch
 133   //
 134
 135   mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
 136   mparam[4]=0; mparam[5]=0; mparam[6]=0;
 137   //
 138   Double_t bparam[6]; // total parameters
 139   Double_t lparam[6]; // local parameters
 140
 141   for (Int_t i=0;i<6;i++) bparam[i]=0;
 142
 143   if (!gGeoManager) {
 144     AliErrorClass("No TGeo\n");
 145     return 0.;
 146   }
 147   //
 148   Double_t length;
 149   Double_t dir[3];
 150   length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
 151                        (end[1]-start[1])*(end[1]-start[1])+
 152                        (end[2]-start[2])*(end[2]-start[2]));
 153   mparam[4]=length;
 154   if (length<TGeoShape::Tolerance()) return 0.0;
 155   Double_t invlen = 1./length;
 156   dir[0] = (end[0]-start[0])*invlen;
 157   dir[1] = (end[1]-start[1])*invlen;
 158   dir[2] = (end[2]-start[2])*invlen;
 159
 160   // Initialize start point and direction
 161   TGeoNode *currentnode = 0;
 162   TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
 163   if (!startnode) {
 164     AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
 165                   start[0],start[1],start[2]));
 166     return 0.0;
 167   }
 168   TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
 169   lparam[0]   = material->GetDensity();
 170   lparam[1]   = material->GetRadLen();
 171   lparam[2]   = material->GetA();
 172   lparam[3]   = material->GetZ();
 173   lparam[4]   = length;
 174   lparam[5]   = lparam[3]/lparam[2];
 175   if (material->IsMixture()) {
 176     TGeoMixture * mixture = (TGeoMixture*)material;
 177     lparam[5] =0;
 178     Double_t sum =0;
 179     for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 180       sum  += mixture->GetWmixt()[iel];
 181       lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 182     }
 183     lparam[5]/=sum;
 184   }
 185
 186   // Locate next boundary within length without computing safety.
 187   // Propagate either with length (if no boundary found) or just cross boundary
 188   gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 189   Double_t step = 0.0; // Step made
 190   Double_t snext = gGeoManager->GetStep();
 191   // If no boundary within proposed length, return current density
 192   if (!gGeoManager->IsOnBoundary()) {
 193     mparam[0] = lparam[0];
 194     mparam[1] = lparam[4]/lparam[1];
 195     mparam[2] = lparam[2];
 196     mparam[3] = lparam[3];
 197     mparam[4] = lparam[4];
 198     return lparam[0];
 199   }
 200   // Try to cross the boundary and see what is next
 201   Int_t nzero = 0;
 202   while (length>TGeoShape::Tolerance()) {
 203     currentnode = gGeoManager->GetCurrentNode();
 204     if (snext<2.*TGeoShape::Tolerance()) nzero++;
 205     else nzero = 0;
 206     if (nzero>3) {
 207       // This means navigation has problems on one boundary
 208       // Try to cross by making a small step
 209       AliErrorClass("Cannot cross boundary\n");
 210       mparam[0] = bparam[0]/step;
 211       mparam[1] = bparam[1];
 212       mparam[2] = bparam[2]/step;
 213       mparam[3] = bparam[3]/step;
 214       mparam[5] = bparam[5]/step;
 215       mparam[4] = step;
 216       mparam[0] = 0.;             // if crash of navigation take mean density 0
 217       mparam[1] = 1000000;        // and infinite rad length
 218       return bparam[0]/step;
 219     }
 220     mparam[6]+=1.;
 221     step += snext;
 222     bparam[1]    += snext/lparam[1];
 223     bparam[2]    += snext*lparam[2];
 224     bparam[3]    += snext*lparam[3];
 225     bparam[5]    += snext*lparam[5];
 226     bparam[0]    += snext*lparam[0];
 227
 228     if (snext>=length) break;
 229     if (!currentnode) break;
 230     length -= snext;
 231     material = currentnode->GetVolume()->GetMedium()->GetMaterial();
 232     lparam[0] = material->GetDensity();
 233     lparam[1]  = material->GetRadLen();
 234     lparam[2]  = material->GetA();
 235     lparam[3]  = material->GetZ();
 236     lparam[5]   = lparam[3]/lparam[2];
 237     if (material->IsMixture()) {
 238       TGeoMixture * mixture = (TGeoMixture*)material;
 239       lparam[5]=0;
 240       Double_t sum =0;
 241       for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 242         sum+= mixture->GetWmixt()[iel];
 243         lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 244       }
 245       lparam[5]/=sum;
 246     }
 247     gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 248     snext = gGeoManager->GetStep();
 249   }
 250   mparam[0] = bparam[0]/step;
 251   mparam[1] = bparam[1];
 252   mparam[2] = bparam[2]/step;
 253   mparam[3] = bparam[3]/step;
 254   mparam[5] = bparam[5]/step;
 255   return bparam[0]/step;
 256 }
 257
 258
 259 Bool_t
 260 AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo,
 261                              Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Double_t sign){
 262   //----------------------------------------------------------------
 263   //
 264   // Propagates the track to the plane X=xk (cm) using the magnetic field map
 265   // and correcting for the crossed material.
 266   //
 267   // mass     - mass used in propagation - used for energy loss correction
 268   // maxStep  - maximal step for propagation
 269   //
 270   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 271   //
 272   //----------------------------------------------------------------
 273   const Double_t kEpsilon = 0.00001;
 274   Double_t xpos     = track->GetX();
 275   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
 276   //
 277   while ( (xToGo-xpos)*dir > kEpsilon){
 278     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 279     Double_t x    = xpos+step;
 280     Double_t xyz0[3],xyz1[3],param[7];
 281     track->GetXYZ(xyz0);   //starting global position
 282
 283     Double_t bz=GetBz(xyz0); // getting the local Bz
 284
 285     if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE;   // no prolongation
 286     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 287
 288     if (TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
 289     if (!track->PropagateTo(x,bz))  return kFALSE;
 290
 291     MeanMaterialBudget(xyz0,xyz1,param);
 292     Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
 293
 294     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 295     if (rotateTo){
 296       if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 297       track->GetXYZ(xyz0);   // global position
 298       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 299       //
 300       Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 301                sa=TMath::Sin(alphan-track->GetAlpha());
 302       Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 303       Double_t sinNew =  sf*ca - cf*sa;
 304       if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 305       if (!track->Rotate(alphan)) return kFALSE;
 306     }
 307     xpos = track->GetX();
 308   }
 309   return kTRUE;
 310 }
 311
 312 Bool_t
 313 AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track,
 314 Double_t xToGo,
 315                                    Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Double_t sign){
 316   //----------------------------------------------------------------
 317   //
 318   // Propagates the track to the plane X=xk (cm)
 319   // taking into account all the three components of the magnetic field
 320   // and correcting for the crossed material.
 321   //
 322   // mass     - mass used in propagation - used for energy loss correction
 323   // maxStep  - maximal step for propagation
 324   //
 325   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 326   //
 327   //----------------------------------------------------------------
 328   const Double_t kEpsilon = 0.00001;
 329   Double_t xpos     = track->GetX();
 330   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
 331   //
 332   while ( (xToGo-xpos)*dir > kEpsilon){
 333     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 334     Double_t x    = xpos+step;
 335     Double_t xyz0[3],xyz1[3],param[7];
 336     track->GetXYZ(xyz0);   //starting global position
 337
 338     Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
 339
 340     if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE;   // no prolongation
 341     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 342
 343     if (TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
 344     if (!track->PropagateToBxByBz(x,b))  return kFALSE;
 345
 346     MeanMaterialBudget(xyz0,xyz1,param);
 347     Double_t xrho=param[0]*param[4]*sign, xx0=param[1];
 348
 349     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 350     if (rotateTo){
 351       if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 352       track->GetXYZ(xyz0);   // global position
 353       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 354       //
 355       Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 356                sa=TMath::Sin(alphan-track->GetAlpha());
 357       Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 358       Double_t sinNew =  sf*ca - cf*sa;
 359       if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 360       if (!track->Rotate(alphan)) return kFALSE;
 361     }
 362     xpos = track->GetX();
 363   }
 364   return kTRUE;
 365 }
 366
 367 Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track,
 368                                            Double_t mass, Double_t step,
 369                                      const AliExternalTrackParam *backup) {
 370   //
 371   // This function brings the "track" with particle "mass" [GeV]
 372   // to the same local coord. system and the same reference plane as
 373   // of the "backup", doing it in "steps" [cm].
 374   // Then, it calculates the 5D predicted Chi2 for these two tracks
 375   //
 376   Double_t chi2=kVeryBig;
 377   Double_t alpha=backup->GetAlpha();
 378   if (!track->Rotate(alpha)) return chi2;
 379
 380   Double_t xb=backup->GetX();
 381   Double_t sign=(xb < track->GetX()) ? 1. : -1.;
 382   if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
 383
 384   chi2=track->GetPredictedChi2(backup);
 385
 386   return chi2;
 387 }
 388
 389
 390
 391
 392 Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1,
 393                    Double_t x2,Double_t y2,
 394                    Double_t x3,Double_t y3)
 395 {
 396   //-----------------------------------------------------------------
 397   // Initial approximation of the track curvature
 398   //-----------------------------------------------------------------
 399   x3 -=x1;
 400   x2 -=x1;
 401   y3 -=y1;
 402   y2 -=y1;
 403   //
 404   Double_t det = x3*y2-x2*y3;
 405   if (TMath::Abs(det)<1e-10) {
 406     return 0;
 407   }
 408   //
 409   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
 410   Double_t x0 = x3*0.5-y3*u;
 411   Double_t y0 = y3*0.5+x3*u;
 412   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
 413   if (det>0) c2*=-1;
 414   return c2;
 415 }
 416
 417 Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1,
 418                    Double_t x2,Double_t y2,
 419                    Double_t x3,Double_t y3)
 420 {
 421   //-----------------------------------------------------------------
 422   // Initial approximation of the track snp
 423   //-----------------------------------------------------------------
 424   x3 -=x1;
 425   x2 -=x1;
 426   y3 -=y1;
 427   y2 -=y1;
 428   //
 429   Double_t det = x3*y2-x2*y3;
 430   if (TMath::Abs(det)<1e-10) {
 431     return 0;
 432   }
 433   //
 434   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
 435   Double_t x0 = x3*0.5-y3*u;
 436   Double_t y0 = y3*0.5+x3*u;
 437   Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0);
 438   x0*=c2;
 439   x0=TMath::Abs(x0);
 440   if (y2*x2<0.) x0*=-1;
 441   return x0;
 442 }
 443
 444 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
 445                    Double_t x2,Double_t y2,
 446                    Double_t z1,Double_t z2, Double_t c)
 447 {
 448   //-----------------------------------------------------------------
 449   // Initial approximation of the tangent of the track dip angle
 450   //-----------------------------------------------------------------
 451   //
 452   x2-=x1;
 453   y2-=y1;
 454   z2-=z1;
 455   Double_t d  =  TMath::Sqrt(x2*x2+y2*y2);  // distance  straight line
 456   if (TMath::Abs(d*c*0.5)>1) return 0;
 457   Double_t   angle2    = TMath::ASin(d*c*0.5);
 458   angle2  = z2*TMath::Abs(c/(angle2*2.));
 459   return angle2;
 460 }
 461
 462
 463 Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
 464                    Double_t x2,Double_t y2,
 465                    Double_t z1,Double_t z2)
 466 {
 467   //-----------------------------------------------------------------
 468   // Initial approximation of the tangent of the track dip angle
 469   //-----------------------------------------------------------------
 470   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 471 }
 472
 473
 474 AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){
 475   //
 476   // Make Seed  - AliExternalTrackParam from input 3 points
 477   // returning seed in local frame of point0
 478   //
 479   Double_t xyz0[3]={0,0,0};
 480   Double_t xyz1[3]={0,0,0};
 481   Double_t xyz2[3]={0,0,0};
 482   Double_t alpha=point0.GetAngle();
 483   Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()};
 484   Double_t bxyz[3]; GetBxByBz(xyz,bxyz);
 485   Double_t bz = bxyz[2];
 486   //
 487   // get points in frame of point 0
 488   //
 489   AliTrackPoint p0r = point0.Rotate(alpha);
 490   AliTrackPoint p1r = point1.Rotate(alpha);
 491   AliTrackPoint p2r = point2.Rotate(alpha);
 492   xyz0[0]=p0r.GetX();
 493   xyz0[1]=p0r.GetY();
 494   xyz0[2]=p0r.GetZ();
 495   xyz1[0]=p1r.GetX();
 496   xyz1[1]=p1r.GetY();
 497   xyz1[2]=p1r.GetZ();
 498   xyz2[0]=p2r.GetX();
 499   xyz2[1]=p2r.GetY();
 500   xyz2[2]=p2r.GetZ();
 501   //
 502   // make covariance estimate
 503   //
 504   Double_t covar[15];
 505   Double_t param[5]={0,0,0,0,0};
 506   for (Int_t m=0; m<15; m++) covar[m]=0;
 507   //
 508   // calculate intitial param
 509   param[0]=xyz0[1];
 510   param[1]=xyz0[2];
 511   param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
 512   param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
 513   param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]);
 514
 515   //covariance matrix - only diagonal elements
 516   //Double_t dist=p0r.GetX()-p2r.GetX();
 517   Double_t deltaP=0;
 518   covar[0]= p0r.GetCov()[3];
 519   covar[2]= p0r.GetCov()[5];
 520   //sigma snp
 521   deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]);
 522   covar[5]+= deltaP*deltaP;
 523   deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]);
 524   covar[5]+= deltaP*deltaP;
 525   deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]);
 526   covar[5]+= deltaP*deltaP;
 527   //sigma tgl
 528   //
 529   deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3];
 530   covar[9]+= deltaP*deltaP;
 531   deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3];
 532   covar[9]+= deltaP*deltaP;
 533   //
 534
 535   deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4];
 536   covar[14]+= deltaP*deltaP;
 537   deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4];
 538   covar[14]+= deltaP*deltaP;
 539   deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4];
 540   covar[14]+= deltaP*deltaP;
 541
 542   covar[14]/=(bz*kB2C)*(bz*kB2C);
 543   param[4]/=(bz*kB2C); // transform to 1/pt
 544   AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar);
 545   if (0) {
 546     // consistency check  -to put warnings here
 547     // small disagrement once Track extrapolation used
 548     // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed
 549     // to check later
 550     Double_t y1,y2,z1,z2;
 551     trackParam->GetYAt(xyz1[0],bz,y1);
 552     trackParam->GetZAt(xyz1[0],bz,z1);
 553     trackParam->GetYAt(xyz2[0],bz,y2);
 554     trackParam->GetZAt(xyz2[0],bz,z2);
 555     if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){
 556       AliWarningClass("Seeding problem y1\n");
 557     }
 558     if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){
 559       AliWarningClass("Seeding problem y2\n");
 560     }
 561     if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){
 562       AliWarningClass("Seeding problem z1\n");
 563     }
 564   }
 565   return trackParam;
 566 }
 567
 568 Double_t  AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){
 569   //
 570   // refit the track  - trackParam using the points in point array
 571   //
 572   const Double_t kMaxSnp=0.99;
 573   if (!trackParam) return 0;
 574   Int_t  npoints=pointArray->GetNPoints();
 575   AliTrackPoint point,point2;
 576   Double_t pointPos[2]={0,0};
 577   Double_t pointCov[3]={0,0,0};
 578   // choose coordinate frame
 579   // in standard way the coordinate frame should be changed point by point
 580   // Some problems with rotation observed
 581   // rotate method of AliExternalTrackParam should be revisited
 582   pointArray->GetPoint(point,0);
 583   pointArray->GetPoint(point2,npoints-1);
 584   Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX());
 585
 586   for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){
 587     pointArray->GetPoint(point,ipoint);
 588     AliTrackPoint pr = point.Rotate(alpha);
 589     trackParam->Rotate(alpha);
 590     Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp);
 591     if(!status){
 592       AliWarningClass("Problem to propagate\n");
 593       break;
 594     }
 595     if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){
 596       AliWarningClass("sin(phi) > kMaxSnp \n");
 597       break;
 598     }
 599     pointPos[0]=pr.GetY();//local y
 600     pointPos[1]=pr.GetZ();//local z
 601     pointCov[0]=pr.GetCov()[3];//simay^2
 602     pointCov[1]=pr.GetCov()[4];//sigmayz
 603     pointCov[2]=pr.GetCov()[5];//sigmaz^2
 604     trackParam->Update(pointPos,pointCov);
 605   }
 606   return 0;
 607 }