X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=STEER%2FAliTrackFitterRieman.cxx;h=3406252e1b7d728718ff3badbd03dc1d7a641208;hb=83e3e2b69c18febd098d851083a4aa2982ebc39b;hp=7fdc505125ab14b7b348592693710813fd179a66;hpb=46ae650f146b19479150f90101be07ef99230470;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliTrackFitterRieman.cxx b/STEER/AliTrackFitterRieman.cxx index 7fdc505125a..3406252e1b7 100644 --- a/STEER/AliTrackFitterRieman.cxx +++ b/STEER/AliTrackFitterRieman.cxx @@ -1,122 +1,211 @@ -#include "TMatrixDSym.h" -#include "TMatrixD.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 + 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"); } //_____________________________________________________________________________ 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)); + fMinPointRadius = rieman.fMinPointRadius; + fMaxPointRadius = rieman.fMaxPointRadius; + 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(UShort_t volId,UShort_t volIdFit, - AliAlignObj::ELayerID layerRangeMin, - AliAlignObj::ELayerID layerRangeMax) +Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, + AliGeomManager::ELayerID layerRangeMin, + AliGeomManager::ELayerID layerRangeMax) { // Fit the track points. The method takes as an input - // the id (volid) of the volume to be skipped from fitting. - // The following two parameters are used to define the + // the set of id's (volids) of the volumes in which + // one wants to calculate the residuals. + // The following parameters are used to define the // range of volumes to be used in the fitting // As a result two AliTrackPointArray's obects are filled. // The first one contains the space points with - // volume id = volid. The second array of points represents - // the track extrapolation corresponding to the space points + // volume id's from volids list. The second array of points represents + // the track extrapolations corresponding to the space points // in the first array. The two arrays can be used to find - // the residuals in the volid and consequently construct a + // the residuals in the volids and consequently construct a // chi2 function to be minimized during the alignment // procedures. For the moment the track extrapolation is taken - // as follows: in XY plane - at the CDA between track circle - // and the space point; in Z - the track extrapolation on the Z - // plane defined by the space point. + // 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,0); - fAlpha = TMath::ATan2(p.GetY(),p.GetX()); +// fPoints->GetPoint(p,0); +// fAlpha = TMath::ATan2(p.GetY(),p.GetX()); Int_t npVolId = 0; fNUsed = 0; @@ -125,32 +214,43 @@ Bool_t AliTrackFitterRieman::Fit(UShort_t volId,UShort_t volIdFit, { fPoints->GetPoint(p,ipoint); UShort_t iVolId = p.GetVolumeID(); - if (iVolId == volId) { + if (FindVolId(volIds,iVolId)) { pindex[npVolId] = ipoint; npVolId++; } - if (volIdFit != 0) { - if (iVolId != volIdFit) continue; + if (volIdsFit != 0x0) { + 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())>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- -// 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 - // + // add point to rieman fitter // - // 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 (0) { - 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 (0) { - TMatrixDSym smatrixz(3); - smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2]; - smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3]; - smatrixz(2,2) = fSumXZ[4]; - smatrixz.Invert(); - TMatrixD sumsxz(1,3); - if (smatrixz.IsValid()){ - sumsxz(0,0) = fSumXZ[5]; - sumsxz(0,1) = fSumXZ[6]; - sumsxz(0,2) = fSumXZ[7]; - TMatrixD res = sumsxz*smatrixz; - fParams[3] = res(0,0); - fParams[4] = res(0,1); - fParams[5] = res(0,2); - for (Int_t i=0;i<3;i++) - for (Int_t j=0;j<=i;j++){ - (*fCov)(i+2,j+2)=smatrixz(i,j); - } - TMatrixD tmp = res*sumsxz.T(); - fChi2 += fSumZZ - tmp(0,0); - fNdf += fNUsed - 3; - conv++; - } - } - else { - Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); - fSumXZ[1] += fSumXZ[5]*Rm1*Rm1; - fSumXZ[2] += fSumXZ[6]*Rm1*Rm1 + fSumXZ[8]*Rm1*Rm1*Rm1*Rm1; - fSumXZ[4] += fSumXZ[7]*Rm1*Rm1; - - TMatrixDSym smatrixz(2); - smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2]; - smatrixz.Invert(); - TMatrixD sumsxz(1,2); - if (smatrixz.IsValid()){ - sumsxz(0,0) = fSumXZ[3]; - sumsxz(0,1) = fSumXZ[4]; - TMatrixD res = sumsxz*smatrixz; - fParams[3] = res(0,0); - fParams[4] = res(0,1); - fParams[5] = 0; - for (Int_t i=0;i<2;i++) - for (Int_t j=0;j<=i;j++){ - (*fCov)(i+3,j+3)=smatrixz(i,j); - } - TMatrixD tmp = res*sumsxz.T(); - fChi2 += fSumZZ - tmp(0,0); - fNdf += fNUsed - 2; - conv++; + 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; } - - // (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){ - 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){ - 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; -} + 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::GetDZat(Double_t 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; + return kTRUE; } -Double_t AliTrackFitterRieman::GetC(){ - return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); -} - -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 @@ -486,27 +430,143 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c // fParam[1] = -x0/y0 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 if (fParams[0] == 0) return kFALSE; + // Track parameters in the global coordinate system 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]; - Double_t x = p.GetX(); - Double_t y = p.GetY(); - Double_t dR = TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)); - Double_t xprime = TMath::Abs(R)*(x-x0)/dR + x0; - Double_t yprime = TMath::Abs(R)*(y-y0)/dR + y0; + // Define space-point refence plane + Double_t alphap = p.GetAngle(); + Double_t sinp = TMath::Sin(alphap); + Double_t cosp = TMath::Cos(alphap); + Double_t x = p.GetX()*cosp + p.GetY()*sinp; + 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)); + return kFALSE; + } + 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; + if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2; + + // Back to the global coordinate system + Double_t xsecond = x*cosp - yprime*sinp; + Double_t ysecond = yprime*cosp + x*sinp; + + // Now Z coordinate and track angles + Double_t x2 = xsecond*cos + ysecond*sin; + Double_t zsecond = GetZat(x2); + Double_t dydx = GetDYat(x2); + Double_t dzdx = GetDZat(x2); + + // Fill the cov matrix of the track extrapolation point + Double_t cov[6] = {0,0,0,0,0,0}; + Double_t sigmax = 100*100.; + cov[0] = sigmax; cov[1] = sigmax*dydx; cov[2] = sigmax*dzdx; + 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+ + cov[3]*sin*sin; + newcov[1] = cov[1]*(cos*cos-sin*sin)- + (cov[3]-cov[0])*sin*cos; + newcov[2] = cov[2]*cos- + cov[4]*sin; + newcov[3] = cov[0]*sin*sin+ + 2*cov[1]*sin*cos+ + cov[3]*cos*cos; + newcov[4] = cov[4]*cos+ + cov[2]*sin; + 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; +} - p2.SetXYZ(xprime,yprime,zprime); +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; +} - return kTRUE; +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; }