X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliTrackFitterRieman.cxx;h=f4e58101ca2f15b1b9eb2802f4f94400ef07e421;hb=313af949bf629682ac449f0b9b760a57808961d1;hp=3b87a1d20ec97f5336299caa338c120175d5a5dc;hpb=24d4520daf9fa97d47602a8729e107dd13e047d7;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliTrackFitterRieman.cxx b/STEER/AliTrackFitterRieman.cxx index 3b87a1d20ec..f4e58101ca2 100644 --- a/STEER/AliTrackFitterRieman.cxx +++ b/STEER/AliTrackFitterRieman.cxx @@ -1,98 +1,142 @@ -#include "TMatrixDSym.h" -#include "TMatrixD.h" -#include "AliTrackFitterRieman.h" +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +/////////////////////////////////////////////////////////////////////////////// +// +// Class to the track points on the Riemann sphere. Inputs are +// the set of id's (volids) of the volumes in which residuals are +// calculated to construct a chi2 function to be minimized during +// the alignment procedures. For the moment the track extrapolation is +// taken 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. +// +// Internal usage of AliRieman class for minimization +// +////////////////////////////////////////////////////////////////////////////// + +#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 + 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 + 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))), + 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"); } //_____________________________________________________________________________ AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman) { - // assignment operator + // + // Assignment operator // 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)); + fDebugStream = 0; + if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root"); return *this; } -AliTrackFitterRieman::~AliTrackFitterRieman() -{ - // destructor + +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 @@ -110,11 +154,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); @@ -141,38 +214,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- - AliAlignObj::kFirstLayer))) 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())>500 || TMath::Abs(plocal.GetX())<2 || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){ + printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<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- -// AliAlignObj::kFirstLayer))) 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 - // - //------------------------------------------------------ - // 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 + // add point to rieman fitter // - 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++; - } + + +Bool_t AliTrackFitterRieman::Update(){ // - // 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;iUpdate(); + 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; } - - // (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, 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){ - 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 { - 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 { - 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 { - 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(){ - 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){ - 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 { + // // Get the closest to a given spacepoint track trajectory point // Look for details in the description of the Fit() method - + // if (!fConv) return kFALSE; // First X and Y coordinates @@ -536,7 +426,7 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c Double_t x0 = -fParams[1]/fParams[0]*cos - 1./fParams[0]*sin; Double_t y0 = 1./fParams[0]*cos - fParams[1]/fParams[0]*sin; if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE; - Double_t R = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/ + Double_t r = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/ fParams[0]; // Define space-point refence plane @@ -547,11 +437,11 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c Double_t y = p.GetY()*cosp - p.GetX()*sinp; Double_t x0p= x0*cosp + y0*sinp; Double_t y0p= y0*cosp - x0*sinp; - if ((R*R - (x-x0p)*(x-x0p))<0) { - 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)); + if ((r*r - (x-x0p)*(x-x0p))<0) { + 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)); return kFALSE; } - Double_t temp = TMath::Sqrt(R*R - (x-x0p)*(x-x0p)); + Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p)); Double_t y1 = y0p + temp; Double_t y2 = y0p - temp; Double_t yprime = y1; @@ -574,6 +464,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+ @@ -590,6 +486,73 @@ 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); +}