STEER/STEER/AliTrackFitterRieman.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$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 // Class to the track points on the Riemann sphere. Inputs are
  21 // the set of id's (volids) of the volumes in which residuals are
  22 // calculated to construct a chi2 function to be minimized during
  23 // the alignment procedures. For the moment the track extrapolation is
  24 // taken at the space-point reference plane. The reference plane is
  25 // found using the covariance matrix of the point
  26 // (assuming sigma(x)=0 at the reference coordinate system.
  27 //
  28 // Internal usage of AliRieman class for minimization
  29 //
  30 //////////////////////////////////////////////////////////////////////////////
  31
  32 #include <TArrayI.h>
  33 #include <TLinearFitter.h>
  34 #include <TMath.h>
  35 #include <TMatrixD.h>
  36 #include <TMatrixDSym.h>
  37 #include <TRandom.h>
  38 #include <TTreeStream.h>
  39
  40 #include "AliLog.h"
  41 #include "AliLog.h"
  42 #include "AliRieman.h"
  43 #include "AliTrackFitterRieman.h"
  44
  45 ClassImp(AliTrackFitterRieman)
  46
  47
  48 AliTrackFitterRieman::AliTrackFitterRieman():
  49   AliTrackFitter(),
  50   fBCorrection(kFALSE),
  51   fAlpha(0.),
  52   fNUsed(0),
  53   fConv(kFALSE),
  54   fMaxDelta(3),
  55   fRieman(new AliRieman(10000)),  // allocate rieman
  56   fMinPointRadius(2.),
  57   fMaxPointRadius(500.),
  58   fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
  59 {
  60   //
  61   // default constructor
  62   //
  63   fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
  64   fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
  65 }
  66
  67
  68 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
  69   AliTrackFitter(array,owner),
  70   fBCorrection(kFALSE),
  71   fAlpha(0.),
  72   fNUsed(0),
  73   fConv(kFALSE),
  74   fMaxDelta(3),
  75   fRieman(new AliRieman(10000)),  //allocate rieman
  76   fMinPointRadius(2.),
  77   fMaxPointRadius(500.),
  78   fDebugStream(0)
  79 {
  80   //
  81   // Constructor
  82   //
  83   fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
  84   fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
  85   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  86 }
  87
  88 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
  89   AliTrackFitter(rieman),
  90   fBCorrection(rieman.fBCorrection),
  91   fAlpha(rieman.fAlpha),
  92   fNUsed(rieman.fNUsed),
  93   fConv(rieman.fConv),
  94   fMaxDelta(rieman.fMaxDelta),
  95   fRieman(new AliRieman(*(rieman.fRieman))),
  96   fMinPointRadius(rieman.fMinPointRadius),
  97   fMaxPointRadius(rieman.fMaxPointRadius),
  98   fDebugStream(0)
  99 {
 100   //
 101   // copy constructor
 102   //
 103   fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
 104   fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
 105   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
 106 }
 107
 108 //_____________________________________________________________________________
 109 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
 110 {
 111   //
 112   // Assignment operator
 113   //
 114   if(this==&rieman) return *this;
 115   ((AliTrackFitter *)this)->operator=(rieman);
 116
 117   fBCorrection = rieman.fBCorrection;
 118   fAlpha  = rieman.fAlpha;
 119   fNUsed  = rieman.fNUsed;
 120   fConv   = rieman.fConv;
 121   fMaxDelta = rieman.fMaxDelta;
 122   fRieman = new AliRieman(*(rieman.fRieman));
 123   fMinPointRadius = rieman.fMinPointRadius;
 124   fMaxPointRadius = rieman.fMaxPointRadius;
 125   fDebugStream = 0;
 126   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
 127   return *this;
 128 }
 129
 130
 131 AliTrackFitterRieman::~AliTrackFitterRieman(){
 132   //
 133   //
 134   //
 135   delete fRieman;
 136   delete fDebugStream;
 137 }
 138
 139 void AliTrackFitterRieman::Reset()
 140 {
 141   // Reset the track parameters and
 142   // rieman sums
 143   //
 144   AliTrackFitter::Reset();
 145   fRieman->Reset();
 146   fAlpha = 0.;
 147   fNUsed = 0;
 148   fConv =kFALSE;
 149 }
 150
 151 Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
 152                                  AliGeomManager::ELayerID layerRangeMin,
 153                          AliGeomManager::ELayerID layerRangeMax)
 154 {
 155   // Fit the track points. The method takes as an input
 156   // the set of id's (volids) of the volumes in which
 157   // one wants to calculate the residuals.
 158   // The following parameters are used to define the
 159   // range of volumes to be used in the fitting
 160   // As a result two AliTrackPointArray's obects are filled.
 161   // The first one contains the space points with
 162   // volume id's from volids list. The second array of points represents
 163   // the track extrapolations corresponding to the space points
 164   // in the first array. The two arrays can be used to find
 165   // the residuals in the volids and consequently construct a
 166   // chi2 function to be minimized during the alignment
 167   // procedures. For the moment the track extrapolation is taken
 168   // at the space-point reference plane. The reference plane is
 169   // found using the covariance matrix of the point
 170   // (assuming sigma(x)=0 at the reference coordinate system.
 171   Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
 172
 173   //  Float_t debugRatio = 1./(1.+debugLevel);
 174   Float_t debugRatio = debugLevel? 1.0/debugLevel : 1.0;
 175
 176   if (!fPoints) {
 177     AliError("Track points array not available! Exiting...");
 178     return kFALSE;
 179   }
 180   Int_t npoints = fPoints->GetNPoints();
 181   if ( npoints<fMinNPoints) return kFALSE;
 182   //
 183   // fast count points
 184   if (volIdsFit != 0x0) {
 185     Int_t countFit = 0;
 186     Int_t countPoint = 0;
 187     for (Int_t ipoint = 0; ipoint < npoints; ipoint++) {
 188       if (FindVolId(volIds,fPoints->GetVolumeID()[ipoint]))
 189         countPoint++;
 190       if (volIdsFit != 0x0) {
 191         if (FindVolId(volIdsFit,fPoints->GetVolumeID()[ipoint]))
 192           countFit++;
 193       }
 194     }
 195     if (countPoint==0) return kFALSE;
 196     if ((countFit<fMinNPoints) && (volIdsFit != 0x0)) return kFALSE;
 197   }
 198   //
 199
 200   Reset();
 201
 202   if (fPoints && volIdsFit && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()<debugRatio){
 203     Int_t nVol    = volIds->GetSize();
 204     Int_t nVolFit = volIdsFit->GetSize();
 205     Int_t volId  = volIds->At(0);
 206     (*fDebugStream)<<"PInput"<<
 207       "NPoints="<<npoints<<    // number of points
 208       "VolId="<<volId<<        // first vol ID
 209       "NVol="<<nVol<<          // number of volumes
 210       "NvolFit="<<nVolFit<<    // number of volumes to fit
 211       "fPoints.="<<fPoints<<   // input points
 212       "\n";
 213   }
 214
 215   Bool_t isAlphaCalc = kFALSE;
 216   AliTrackPoint p,plocal;
 217 //   fPoints->GetPoint(p,0);
 218 //   fAlpha = TMath::ATan2(p.GetY(),p.GetX());
 219
 220   Int_t npVolId = 0;
 221   fNUsed = 0;
 222   Int_t *pindex = new Int_t[npoints];
 223   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 224     {
 225       fPoints->GetPoint(p,ipoint);
 226       UShort_t iVolId = p.GetVolumeID();
 227       if (FindVolId(volIds,iVolId)) {
 228         pindex[npVolId] = ipoint;
 229         npVolId++;
 230       }
 231       if (volIdsFit != 0x0) {
 232         if (!FindVolId(volIdsFit,iVolId)) continue;
 233       }
 234       else {
 235         if (iVolId < AliGeomManager::LayerToVolUID(layerRangeMin,0) ||
 236             iVolId > AliGeomManager::LayerToVolUID(layerRangeMax,
 237                                                 AliGeomManager::LayerSize(layerRangeMax))) continue;
 238       }
 239       if (!isAlphaCalc) {
 240         fAlpha = p.GetAngle();
 241         isAlphaCalc = kTRUE;
 242       }
 243       plocal = p.Rotate(fAlpha);
 244       if (TMath::Abs(plocal.GetX())>fMaxPointRadius || TMath::Abs(plocal.GetX())<fMinPointRadius || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
 245         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
 246         p.Dump();
 247         plocal.Dump();
 248         printf("Problematic point\n");
 249         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
 250       }else{
 251         AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
 252                  TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
 253       }
 254       // fNUsed++;  AddPoint should be responsible
 255     }
 256
 257   if (npVolId == 0 || fNUsed < fMinNPoints) {
 258     delete [] pindex;
 259     return kFALSE;
 260   }
 261
 262   Update();
 263
 264   if (!fConv) {
 265     delete [] pindex;
 266     return kFALSE;
 267   }
 268
 269   if ((fParams[0] == 0) ||
 270       ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
 271     delete [] pindex;
 272     return kFALSE;
 273   }
 274
 275
 276   if (fNUsed < fMinNPoints) {
 277     delete [] pindex;
 278     return kFALSE;
 279   }
 280
 281   fPVolId = new AliTrackPointArray(npVolId);
 282   fPTrack = new AliTrackPointArray(npVolId);
 283   AliTrackPoint p2;
 284   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
 285     {
 286       Int_t index = pindex[ipoint];
 287       fPoints->GetPoint(p,index);
 288       if (GetPCA(p,p2) && (
 289           TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta &&
 290           TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
 291           TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
 292           )) {
 293         Float_t xyz[3],xyz2[3];
 294         p.GetXYZ(xyz); p2.GetXYZ(xyz2);
 295         //      printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);
 296         fPVolId->AddPoint(ipoint,&p);
 297         fPTrack->AddPoint(ipoint,&p2);
 298       }else{
 299         // what should be default bahavior -
 300         delete [] pindex;
 301         delete fPVolId;
 302         delete fPTrack;
 303         fPVolId =0;
 304         fPTrack =0;
 305         return kFALSE;
 306       }
 307     }
 308
 309   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
 310     //
 311     // Debug Info
 312     //
 313     AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
 314     AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
 315     AliTrackPointArray *lPTrackE = new AliTrackPointArray(npVolId);
 316     AliRieman * residual = fRieman->MakeResiduals();
 317     for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
 318       AliTrackPoint p0, p0local;
 319       AliTrackPoint pFit, pFitlocal, pFitLocalE;
 320       fPVolId->GetPoint(p0,ipoint);
 321       Float_t lAngle = p0.GetAngle();
 322       p0local= p0.MasterToLocal();
 323       fPTrack->GetPoint(pFit,ipoint);
 324       pFitlocal= pFit.Rotate(lAngle);
 325       //
 326       Float_t xyz[3], cov[6];
 327       xyz[0] = pFitlocal.GetX();
 328       xyz[1] = pFitlocal.GetY();
 329       xyz[2] = pFitlocal.GetZ();
 330       for (Int_t icov=0; icov<6; icov++) cov[icov]=0;
 331       cov[3] = GetErrY2at(xyz[0]);
 332       cov[5] = GetErrZ2at(xyz[0]);
 333       pFitLocalE.SetXYZ(xyz,cov);
 334       //
 335       lPVolId->AddPoint(ipoint,&p0local);
 336       lPTrack->AddPoint(ipoint,&pFitlocal);
 337       lPTrackE->AddPoint(ipoint,&pFitLocalE);
 338     }
 339     //
 340     // debug info
 341     //
 342     Int_t nVol    = volIds->GetSize();
 343     Int_t nVolFit = volIdsFit->GetSize();
 344     Int_t volId   = volIds->At(0);
 345     Int_t modId   =0;
 346     Int_t layer   =  AliGeomManager::VolUIDToLayer(volId,modId);
 347     Int_t volIdFit   = volIdsFit->At(0);
 348     Int_t modIdFit   =0;
 349     Int_t layerFit   =  AliGeomManager::VolUIDToLayer(volIdFit,modIdFit);
 350
 351     (*fDebugStream)<<"Fit"<<
 352       "VolId="<<volId<<        // volume ID
 353       "Layer="<<layer<<        // layer ID
 354       "Module="<<modId<<       // module ID
 355       "LayerFit="<<layerFit<<        // layer ID fit
 356       "ModuleFit="<<modIdFit<<       // module ID fit
 357       "NVol="<<nVol<<          // number of volumes
 358       "NvolFit="<<nVolFit<<    // number of volumes to fit
 359       "Points0.="<<fPVolId<<    // original points
 360       "Points1.="<<fPTrack<<    // fitted points
 361       "LPoints0.="<<lPVolId<<   // original points - local frame
 362       "LPoints1.="<<lPTrack<<   // fitted points   - local frame
 363       "LPointsE.="<<lPTrackE<<   // fitted points with ext error - local frame
 364       "Rieman.="<<this<<        // original rieman fit
 365       "Res.="<<residual<<       // residuals of rieman fit
 366       "\n";
 367     delete lPVolId;
 368     delete lPTrack;
 369     delete residual;
 370   }
 371
 372   delete [] pindex;
 373   return kTRUE;
 374 }
 375
 376
 377 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
 378 {
 379   //
 380   // add point to rieman fitter
 381   //
 382   fRieman->AddPoint(x,y,z,sy,sz);
 383   fNUsed = fRieman->GetN();
 384 }
 385
 386
 387
 388 Bool_t AliTrackFitterRieman::Update(){
 389   //
 390   //
 391   //
 392   fRieman->Update();
 393   fConv = kFALSE;
 394   if (fRieman->IsValid()){
 395     for (Int_t ipar=0; ipar<6; ipar++){
 396       fParams[ipar] = fRieman->GetParam()[ipar];
 397     }
 398     fChi2  =  fRieman->GetChi2();
 399     fNdf   =  fRieman->GetN()- 2;
 400     fNUsed =  fRieman->GetN();
 401     fConv = kTRUE;
 402   }
 403   //
 404   //
 405   //
 406   TLinearFitter fitY(3,"pol2");
 407   TLinearFitter fitZ(3,"pol2");
 408   for (Int_t ip=0; ip<fRieman->GetN();ip++){
 409     Double_t x = fRieman->GetX()[ip];
 410     fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1);
 411     fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1);
 412   }
 413   fitY.Eval();
 414   fitZ.Eval();
 415   for (Int_t iparam=0; iparam<3; iparam++){
 416     fCorrY[iparam]=fitY.GetParameter(iparam);
 417     fCorrZ[iparam]=fitZ.GetParameter(iparam);
 418   }
 419   fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3);
 420   fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3);
 421
 422   return kTRUE;
 423 }
 424
 425
 426
 427 //_____________________________________________________________________________
 428 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
 429 {
 430   //
 431   // Get the closest to a given spacepoint track trajectory point
 432   // Look for details in the description of the Fit() method
 433   //
 434   if (!fConv) return kFALSE;
 435
 436   // First X and Y coordinates
 437   Double_t sin = TMath::Sin(fAlpha);
 438   Double_t cos = TMath::Cos(fAlpha);
 439   //  fParam[0]  = 1/y0
 440   //  fParam[1]  = -x0/y0
 441   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 442   if (fParams[0] == 0) return kFALSE;
 443   // Track parameters in the global coordinate system
 444   Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
 445   Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
 446   if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
 447   Double_t r  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
 448                 fParams[0];
 449
 450   // Define space-point refence plane
 451   Double_t alphap = p.GetAngle();
 452   Double_t sinp = TMath::Sin(alphap);
 453   Double_t cosp = TMath::Cos(alphap);
 454   Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
 455   Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
 456   Double_t x0p= x0*cosp + y0*sinp;
 457   Double_t y0p= y0*cosp - x0*sinp;
 458   if ((r*r - (x-x0p)*(x-x0p))<0) {
 459     AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",r,x0p,x,alphap));
 460     return kFALSE;
 461   }
 462   Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
 463   Double_t y1 = y0p + temp;
 464   Double_t y2 = y0p - temp;
 465   Double_t yprime = y1;
 466   if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
 467
 468   // Back to the global coordinate system
 469   Double_t xsecond = x*cosp - yprime*sinp;
 470   Double_t ysecond = yprime*cosp + x*sinp;
 471
 472   // Now Z coordinate and track angles
 473   Double_t x2 = xsecond*cos + ysecond*sin;
 474   Double_t zsecond = GetZat(x2);
 475   Double_t dydx = GetDYat(x2);
 476   Double_t dzdx = GetDZat(x2);
 477
 478   // Fill the cov matrix of the track extrapolation point
 479   Double_t cov[6] = {0,0,0,0,0,0};
 480   Double_t sigmax = 100*100.;
 481   cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
 482   cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
 483   cov[5] = sigmax*dzdx*dzdx;
 484
 485   Double_t sigmay2 = GetErrY2at(x2);
 486   Double_t sigmaz2 = GetErrZ2at(x2);
 487   cov[3] += sigmay2;
 488   cov[5] += sigmaz2;
 489
 490
 491   Float_t  newcov[6];
 492   newcov[0] = cov[0]*cos*cos-
 493             2*cov[1]*sin*cos+
 494               cov[3]*sin*sin;
 495   newcov[1] = cov[1]*(cos*cos-sin*sin)-
 496              (cov[3]-cov[0])*sin*cos;
 497   newcov[2] = cov[2]*cos-
 498               cov[4]*sin;
 499   newcov[3] = cov[0]*sin*sin+
 500             2*cov[1]*sin*cos+
 501               cov[3]*cos*cos;
 502   newcov[4] = cov[4]*cos+
 503               cov[2]*sin;
 504   newcov[5] = cov[5];
 505
 506   p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
 507   Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
 508   Float_t debugRatio = 1./(1.+debugLevel);
 509   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
 510     AliTrackPoint lp0(p);
 511     AliTrackPoint lp2(p2);
 512     AliTrackPoint localp0(p);
 513     AliTrackPoint localp2(p2);
 514     Float_t lAngle = lp0.GetAngle();
 515     localp0 = localp0.Rotate(lAngle);
 516     localp2 = localp2.Rotate(lAngle);
 517
 518     (*fDebugStream)<<"PCA"<<
 519       "P0.="<<&lp0<<  //global position
 520       "P2.="<<&lp2<<
 521       "LP0.="<<&localp0<<  //local position
 522       "LP2.="<<&localp2<<
 523       "\n";
 524   }
 525   return kTRUE;
 526 }
 527
 528 Double_t AliTrackFitterRieman::GetYat(Double_t x) const  {
 529   //
 530   // get y position at given point
 531   //
 532   Double_t correction=0;
 533   if (fBCorrection){ // systematic effect correction
 534     correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
 535   }
 536   return fRieman->GetYat(x)+correction;
 537 }
 538
 539 Double_t AliTrackFitterRieman::GetZat(Double_t x) const  {
 540   //
 541   // get z position at given point
 542   //
 543   Double_t correction=0;
 544   if (fBCorrection){ // systematic effect correction
 545     correction  =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
 546   }
 547   return fRieman->GetZat(x)+correction;
 548 }
 549
 550 Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const {
 551   //
 552   // get estimate of extrapolation error
 553   //
 554   Double_t error = fRieman->GetErrY(x);
 555   Double_t correction=0;
 556   if (fBCorrection){ // everestimate error due systematic effect
 557     error      *=fCorrY[3];
 558     correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
 559     correction *=correction;
 560   }
 561   return TMath::Sqrt(error+correction);
 562 }
 563
 564 Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const {
 565   //
 566   // get estimate of extrapolation error
 567   //
 568   Double_t error = fRieman->GetErrZ(x)*fCorrZ[3];
 569   Double_t correction=0;
 570   if (fBCorrection){
 571      error    *=  fCorrZ[3];
 572     correction =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
 573     correction*=  correction;
 574   }
 575   return TMath::Sqrt(error+correction);
 576 }
 577
 578 void AliTrackFitterRieman::SetParam(Int_t i, Double_t par) {
 579   if (i<0 || i>5) return;
 580   fParams[i]=par;
 581   fRieman->GetParam()[i]=par;
 582 }