X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliTrackFitterRieman.cxx;h=3406252e1b7d728718ff3badbd03dc1d7a641208;hb=25e05e8a7c9beecd23fb055239709a2bde4b64ec;hp=112f0be92a954e44e56ae4dcb76604b47646f5ca;hpb=6b6cba33c24c613006a073283a9fbbfb58f79032;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliTrackFitterRieman.cxx b/STEER/AliTrackFitterRieman.cxx index 112f0be92a9..3406252e1b7 100644 --- a/STEER/AliTrackFitterRieman.cxx +++ b/STEER/AliTrackFitterRieman.cxx @@ -25,59 +25,78 @@ // found using the covariance matrix of the point // (assuming sigma(x)=0 at the reference coordinate system. // +// Internal usage of AliRieman class for minimization +// ////////////////////////////////////////////////////////////////////////////// -#include "TMatrixDSym.h" -#include "TMatrixD.h" -#include "AliTrackFitterRieman.h" +#include +#include +#include +#include +#include +#include +#include + #include "AliLog.h" +#include "AliLog.h" +#include "AliRieman.h" +#include "AliTrackFitterRieman.h" ClassImp(AliTrackFitterRieman) + AliTrackFitterRieman::AliTrackFitterRieman(): - AliTrackFitter() + AliTrackFitter(), + fBCorrection(kFALSE), + fAlpha(0.), + fNUsed(0), + fConv(kFALSE), + fMaxDelta(3), + fRieman(new AliRieman(10000)), // allocate rieman + fMinPointRadius(2.), + fMaxPointRadius(500.), + fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root")) { // // default constructor // - fAlpha = 0.; - for (Int_t i=0;i<9;i++) fSumXY[i] = 0; - fSumYY = 0; - for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; - fSumZZ = 0; - fNUsed = 0; - fConv = kFALSE; } AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner): - AliTrackFitter(array,owner) + AliTrackFitter(array,owner), + fBCorrection(kFALSE), + fAlpha(0.), + fNUsed(0), + fConv(kFALSE), + fMaxDelta(3), + fRieman(new AliRieman(10000)), //allocate rieman + fMinPointRadius(2.), + fMaxPointRadius(500.), + fDebugStream(0) { // // Constructor // - fAlpha = 0.; - for (Int_t i=0;i<9;i++) fSumXY[i] = 0; - fSumYY = 0; - for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; - fSumZZ = 0; - fNUsed = 0; - fConv = kFALSE; + if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root"); } AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman): - AliTrackFitter(rieman) + AliTrackFitter(rieman), + fBCorrection(rieman.fBCorrection), + fAlpha(rieman.fAlpha), + fNUsed(rieman.fNUsed), + fConv(rieman.fConv), + fMaxDelta(rieman.fMaxDelta), + fRieman(new AliRieman(*(rieman.fRieman))), + fMinPointRadius(rieman.fMinPointRadius), + fMaxPointRadius(rieman.fMaxPointRadius), + fDebugStream(0) { // // copy constructor // - fAlpha = rieman.fAlpha; - for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i]; - fSumYY = rieman.fSumYY; - for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i]; - fSumZZ = rieman.fSumZZ; - fNUsed = rieman.fNUsed; - fConv = rieman.fConv; + if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root"); } //_____________________________________________________________________________ @@ -89,37 +108,43 @@ AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRiema if(this==&rieman) return *this; ((AliTrackFitter *)this)->operator=(rieman); - fAlpha = rieman.fAlpha; - for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i]; - fSumYY = rieman.fSumYY; - for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i]; - fSumZZ = rieman.fSumZZ; - fNUsed = rieman.fNUsed; - fConv = rieman.fConv; - + fBCorrection = rieman.fBCorrection; + fAlpha = rieman.fAlpha; + fNUsed = rieman.fNUsed; + fConv = rieman.fConv; + fMaxDelta = rieman.fMaxDelta; + fRieman = new AliRieman(*(rieman.fRieman)); + fMinPointRadius = rieman.fMinPointRadius; + fMaxPointRadius = rieman.fMaxPointRadius; + fDebugStream = 0; + if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root"); return *this; } -//_____________________________________________________________________________ + +AliTrackFitterRieman::~AliTrackFitterRieman(){ + // + // + // + delete fRieman; + delete fDebugStream; +} + void AliTrackFitterRieman::Reset() { - // // Reset the track parameters and // rieman sums // AliTrackFitter::Reset(); + fRieman->Reset(); fAlpha = 0.; - for (Int_t i=0;i<9;i++) fSumXY[i] = 0; - fSumYY = 0; - for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; - fSumZZ = 0; fNUsed = 0; fConv =kFALSE; } Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, - AliAlignObj::ELayerID layerRangeMin, - AliAlignObj::ELayerID layerRangeMax) + AliGeomManager::ELayerID layerRangeMin, + AliGeomManager::ELayerID layerRangeMax) { // Fit the track points. The method takes as an input // the set of id's (volids) of the volumes in which @@ -137,11 +162,45 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, // at the space-point reference plane. The reference plane is // found using the covariance matrix of the point // (assuming sigma(x)=0 at the reference coordinate system. + Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman"); + + // Float_t debugRatio = 1./(1.+debugLevel); + Float_t debugRatio = debugLevel? 1.0/debugLevel : 1.0; + + Int_t npoints = fPoints->GetNPoints(); + if ( npointsGetVolumeID()[ipoint])) + countPoint++; + if (volIdsFit != 0x0) { + if (FindVolId(volIdsFit,fPoints->GetVolumeID()[ipoint])) + countFit++; + } + } + if (countPoint==0) return kFALSE; + if ((countFitGetNPoints(); - if (npoints < 3) return kFALSE; + if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()GetSize(); + Int_t nVolFit = volIdsFit->GetSize(); + Int_t volId = volIds->At(0); + (*fDebugStream)<<"PInput"<< + "NPoints="<GetPoint(p,ipoint); @@ -168,37 +222,34 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, if (!FindVolId(volIdsFit,iVolId)) continue; } else { - if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) || - iVolId > AliAlignObj::LayerToVolUID(layerRangeMax, - AliAlignObj::LayerSize(layerRangeMax))) continue; + if (iVolId < AliGeomManager::LayerToVolUID(layerRangeMin,0) || + iVolId > AliGeomManager::LayerToVolUID(layerRangeMax, + AliGeomManager::LayerSize(layerRangeMax))) continue; } if (!isAlphaCalc) { fAlpha = p.GetAngle(); isAlphaCalc = kTRUE; } plocal = p.Rotate(fAlpha); - AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(), - TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5])); - fNUsed++; + if (TMath::Abs(plocal.GetX())>fMaxPointRadius || TMath::Abs(plocal.GetX())GetPoint(p,index); - if (GetPCA(p,p2)) { + if (GetPCA(p,p2) && ( + TMath::Abs(p.GetX()-p2.GetX())AddPoint(ipoint,&p); fPTrack->AddPoint(ipoint,&p2); + }else{ + // what should be default bahavior - + delete [] pindex; + delete fPVolId; + delete fPTrack; + fPVolId =0; + fPTrack =0; + return kFALSE; } } + + if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()MakeResiduals(); + for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){ + AliTrackPoint p0, p0local; + AliTrackPoint pFit, pFitlocal, pFitLocalE; + fPVolId->GetPoint(p0,ipoint); + Float_t lAngle = p0.GetAngle(); + p0local= p0.MasterToLocal(); + fPTrack->GetPoint(pFit,ipoint); + pFitlocal= pFit.Rotate(lAngle); + // + Float_t xyz[3], cov[6]; + xyz[0] = pFitlocal.GetX(); + xyz[1] = pFitlocal.GetY(); + xyz[2] = pFitlocal.GetZ(); + for (Int_t icov=0; icov<6; icov++) cov[icov]=0; + cov[3] = GetErrY2at(xyz[0]); + cov[5] = GetErrZ2at(xyz[0]); + pFitLocalE.SetXYZ(xyz,cov); + // + lPVolId->AddPoint(ipoint,&p0local); + lPTrack->AddPoint(ipoint,&pFitlocal); + lPTrackE->AddPoint(ipoint,&pFitLocalE); + } + // + // debug info + // + Int_t nVol = volIds->GetSize(); + Int_t nVolFit = volIdsFit->GetSize(); + Int_t volId = volIds->At(0); + Int_t modId =0; + Int_t layer = AliGeomManager::VolUIDToLayer(volId,modId); + Int_t volIdFit = volIdsFit->At(0); + Int_t modIdFit =0; + Int_t layerFit = AliGeomManager::VolUIDToLayer(volIdFit,modIdFit); + + (*fDebugStream)<<"Fit"<< + "VolId="<GetPoint(p,ipoint); -// UShort_t iVolId = p.GetVolumeID(); -// if (volIdFit != 0) { -// if (iVolId != volIdFit) continue; -// } -// else { -// if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) || -// iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,AliAlignObj::LayerSize(layerRangeMax))) continue; -// } -// plocal = p.Rotate(fAlpha); -// Float_t delta = (fParams[0]*(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY())+ -// 2.*plocal.GetX()*fParams[1]+ -// fParams[2]- -// 2.*plocal.GetY())/ -// (2.*TMath::Sqrt(plocal.GetCov()[3])); -// // Float_t delta2 = (fParams[3]+ -// // plocal.GetX()*fParams[4]+ -// // plocal.GetX()*plocal.GetX()*fParams[5]- -// // plocal.GetZ())/ -// // (TMath::Sqrt(plocal.GetCov()[5])); -// Double_t r = TMath::Sqrt(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY()); -// Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); -// Float_t delta2 = (fParams[3]+ -// r*fParams[4]+r*r*r*fParams[4]*Rm1*Rm1/24- -// plocal.GetZ())/ -// (TMath::Sqrt(plocal.GetCov()[5])); -// chi2 += delta*delta; -// chi22 += delta2*delta2; -// // printf("pulls %d %d %f %f\n",ipoint,iVolId,delta,delta2); - -// } -// printf("My chi2 = %f + %f = %f\n",chi2,chi22,chi2+chi22); - return kTRUE; } + void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz) { // - // Rieman update + // add point to rieman fitter // - //------------------------------------------------------ - // XY direction - // - // (x-x0)^2+(y-y0)^2-R^2=0 ===> - // - // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> - // - // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2 - // - // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation - // - // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 - // - // final linear equation : a + u*b +t*c - v =0; - // - // Minimization : - // - // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min; - // - // where sigmai is the error of maesurement (a + ui*b +ti*c - vi) - // - // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t - // - fX[fNUsed] = x; fY[fNUsed]=y; fZ[fNUsed]=z; fSy[fNUsed]=sy; fSz[fNUsed]=sz; - // - // XY part - // - Double_t t = x*x+y*y; - if (t<2) return; - t = 1./t; - Double_t u = 2.*x*t; - Double_t v = 2.*y*t; - Double_t error = 2.*sy*t; - error *=error; - Double_t weight = 1./error; - fSumXY[0] +=weight; - fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight; - fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight; - fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight; - fSumYY += v*v*weight; - // - // XZ part - // - if (1) { - weight = 1./(sz*sz); -// fSumXZ[0] +=weight; -// fSumXZ[1] +=weight*x; fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x; -// fSumXZ[5] +=weight*z; fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z; - fSumZZ += z*z*weight; - } - else { - weight = 1./(sz*sz); - fSumXZ[0] +=weight; - Double_t r = TMath::Sqrt(x*x+y*y); - fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z; - // Now the auxulary sums - fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24; - fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24); - fSumZZ += z*z*weight; - } + fRieman->AddPoint(x,y,z,sy,sz); + fNUsed = fRieman->GetN(); } -void AliTrackFitterRieman::Update(){ - // - // Rieman update - // - // - for (Int_t i=0;i<6;i++)fParams[i]=0; - fChi2 = 0; - fNdf = 0; - Int_t conv=0; - // - // XY part - // - TMatrixDSym smatrix(3); - TMatrixD sums(1,3); - // - // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; - // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; - // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; - - smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; - smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; - sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; - smatrix.Invert(); - if (smatrix.IsValid()){ - for (Int_t i=0;i<3;i++) - for (Int_t j=0;j<=i;j++){ - (*fCov)(i,j)=smatrix(i,j); - } - TMatrixD res = sums*smatrix; - fParams[0] = res(0,0); - fParams[1] = res(0,1); - fParams[2] = res(0,2); - TMatrixD tmp = res*sums.T(); - fChi2 += fSumYY - tmp(0,0); - fNdf += fNUsed - 3; - conv++; - } - // - // XZ part - // - if (1) { - Double_t x0 = -fParams[1]/fParams[0]; - Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); - - for (Int_t i=0;i +Bool_t AliTrackFitterRieman::Update(){ // - // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> - // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 + // // - // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation - // - // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 - // final linear equation : a + u*b +t*c - v =0; - // - // fParam[0] = 1/y0 - // fParam[1] = -x0/y0 - // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 - if (conv>1) fConv =kTRUE; - else - fConv=kFALSE; -} - -//_____________________________________________________________________________ -Double_t AliTrackFitterRieman::GetYat(Double_t x) const -{ - // - // Returns value of Y at given X - // - if (!fConv) return 0.; - Double_t res2 = (x*fParams[0]+fParams[1]); - res2*=res2; - res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; - if (res2<0) return 0; - res2 = TMath::Sqrt(res2); - res2 = (1-res2)/fParams[0]; - return res2; -} - -//_____________________________________________________________________________ -Double_t AliTrackFitterRieman::GetDYat(Double_t x) const -{ - // - // Returns value of dY/dX at given X - // - if (!fConv) return 0.; - Double_t x0 = -fParams[1]/fParams[0]; - if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0; - Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); - if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0; - Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0)); - if (fParams[0]<0) res*=-1.; - return res; -} - -//_____________________________________________________________________________ -Double_t AliTrackFitterRieman::GetZat(Double_t x) const -{ - // - // Returns value of Z given X + fRieman->Update(); + fConv = kFALSE; + if (fRieman->IsValid()){ + for (Int_t ipar=0; ipar<6; ipar++){ + fParams[ipar] = fRieman->GetParam()[ipar]; + } + fChi2 = fRieman->GetChi2(); + fNdf = fRieman->GetN()- 2; + fNUsed = fRieman->GetN(); + fConv = kTRUE; + } // - if (!fConv) return 0.; - Double_t x0 = -fParams[1]/fParams[0]; - if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; - Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); - Double_t phi = TMath::ASin((x-x0)*rm1); - Double_t phi0 = TMath::ASin((0.-x0)*rm1); - Double_t dphi = (phi-phi0); - Double_t res = fParams[3]+fParams[4]*dphi/rm1; - return res; -} - -//_____________________________________________________________________________ -Double_t AliTrackFitterRieman::GetDZat(Double_t x) const -{ // - // Returns value of dZ/dX given X // - if (!fConv) return 0.; - Double_t x0 = -fParams[1]/fParams[0]; - if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; - Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); - Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1)); - return res; -} - + TLinearFitter fitY(3,"pol2"); + TLinearFitter fitZ(3,"pol2"); + for (Int_t ip=0; ipGetN();ip++){ + Double_t x = fRieman->GetX()[ip]; + fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1); + fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1); + } + fitY.Eval(); + fitZ.Eval(); + for (Int_t iparam=0; iparam<3; iparam++){ + fCorrY[iparam]=fitY.GetParameter(iparam); + fCorrZ[iparam]=fitZ.GetParameter(iparam); + } + fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3); + fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3); -//_____________________________________________________________________________ -Double_t AliTrackFitterRieman::GetC() const -{ - // - // Returns curvature - // - return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); + return kTRUE; } -//_____________________________________________________________________________ -Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz) const -{ - // - // Returns position given radius - // - if (!fConv) return kFALSE; - Double_t res2 = (r*fParams[0]+fParams[1]); - res2*=res2; - res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; - if (res2<0) return kFALSE; - res2 = TMath::Sqrt(res2); - res2 = (1-res2)/fParams[0]; - - if (!fConv) return kFALSE; - Double_t x0 = -fParams[1]/fParams[0]; - if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; - Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); - Double_t phi = TMath::ASin((r-x0)*rm1); - Double_t phi0 = TMath::ASin((0.-x0)*rm1); - Double_t dphi = (phi-phi0); - Double_t res = fParams[3]+fParams[4]*dphi/rm1; - Double_t sin = TMath::Sin(fAlpha); - Double_t cos = TMath::Cos(fAlpha); - xyz[0] = r*cos - res2*sin; - xyz[1] = res2*cos + r*sin; - xyz[2] = res; - - return kTRUE; -} //_____________________________________________________________________________ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const @@ -629,6 +472,12 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx; cov[5] = sigmax*dzdx*dzdx; + Double_t sigmay2 = GetErrY2at(x2); + Double_t sigmaz2 = GetErrZ2at(x2); + cov[3] += sigmay2; + cov[5] += sigmaz2; + + Float_t newcov[6]; newcov[0] = cov[0]*cos*cos- 2*cov[1]*sin*cos+ @@ -645,6 +494,79 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c newcov[5] = cov[5]; p2.SetXYZ(xsecond,ysecond,zsecond,newcov); - + Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman"); + Float_t debugRatio = 1./(1.+debugLevel); + if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()GetYat(x)+correction; +} + +Double_t AliTrackFitterRieman::GetZat(Double_t x) const { + // + // get z position at given point + // + Double_t correction=0; + if (fBCorrection){ // systematic effect correction + correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x; + } + return fRieman->GetZat(x)+correction; +} + +Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const { + // + // get estimate of extrapolation error + // + Double_t error = fRieman->GetErrY(x); + Double_t correction=0; + if (fBCorrection){ // everestimate error due systematic effect + error *=fCorrY[3]; + correction = fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x; + correction *=correction; + } + return TMath::Sqrt(error+correction); +} + +Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const { + // + // get estimate of extrapolation error + // + Double_t error = fRieman->GetErrZ(x)*fCorrZ[3]; + Double_t correction=0; + if (fBCorrection){ + error *= fCorrZ[3]; + correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x; + correction*= correction; + } + return TMath::Sqrt(error+correction); +} + +void AliTrackFitterRieman::SetParam(Int_t i, Double_t par) { + if (i<0 || i>5) return; + fParams[i]=par; + fRieman->GetParam()[i]=par; +}