From dfce6628827b03c8c76b32843c373692217f7209 Mon Sep 17 00:00:00 2001 From: cvetan Date: Tue, 26 Sep 2006 08:38:20 +0000 Subject: [PATCH] New versions of rieman track fitter and residuals minimization based on root TLinearFitter (M.Ivanov). Some effective C++ mods --- STEER/AliTrackFitterRieman.cxx | 164 ++++++++++++++++++++++++++---- STEER/AliTrackFitterRieman.h | 16 ++- STEER/AliTrackResidualsLinear.cxx | 56 ++++++++-- STEER/AliTrackResidualsLinear.h | 5 +- 4 files changed, 205 insertions(+), 36 deletions(-) diff --git a/STEER/AliTrackFitterRieman.cxx b/STEER/AliTrackFitterRieman.cxx index 583cd064469..7a9a3c52f33 100644 --- a/STEER/AliTrackFitterRieman.cxx +++ b/STEER/AliTrackFitterRieman.cxx @@ -32,17 +32,20 @@ #include "TMatrixDSym.h" #include "TMatrixD.h" #include "TArrayI.h" +#include "TLinearFitter.h" #include "AliTrackFitterRieman.h" #include "AliLog.h" #include "TTreeStream.h" #include "AliRieman.h" #include "AliLog.h" +#include "TRandom.h" ClassImp(AliTrackFitterRieman) AliTrackFitterRieman::AliTrackFitterRieman(): AliTrackFitter(), + fBCorrection(kFALSE), fAlpha(0.), fNUsed(0), fConv(kFALSE), @@ -58,6 +61,7 @@ AliTrackFitterRieman::AliTrackFitterRieman(): AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner): AliTrackFitter(array,owner), + fBCorrection(kFALSE), fAlpha(0.), fNUsed(0), fConv(kFALSE), @@ -73,6 +77,7 @@ AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t own AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman): AliTrackFitter(rieman), + fBCorrection(rieman.fBCorrection), fAlpha(rieman.fAlpha), fNUsed(rieman.fNUsed), fConv(rieman.fConv), @@ -95,6 +100,7 @@ AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRiema if(this==&rieman) return *this; ((AliTrackFitter *)this)->operator=(rieman); + fBCorrection = rieman.fBCorrection; fAlpha = rieman.fAlpha; fNUsed = rieman.fNUsed; fConv = rieman.fConv; @@ -146,10 +152,35 @@ 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); - Reset(); + const Int_t kMinPoints =1; Int_t npoints = fPoints->GetNPoints(); - if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){ + if ( npointsGetSize(); ifit++){ + Int_t volIdFit = volIdsFit->At(ifit); + for (Int_t ipoint = 0; ipoint < npoints; ipoint++) + if (volIdFit==fPoints->GetVolumeID()[ipoint]) countFit++; + } + if (countFitGetSize(); jpoint++){ + Int_t volIdPoint = volIds->At(jpoint); + for (Int_t ipoint = 0; ipoint < npoints; ipoint++) + if (volIdPoint==fPoints->GetVolumeID()[ipoint]) countPoint++; + } + if (countPoint1&& gRandom->Rndm()GetSize(); Int_t nVolFit = volIdsFit->GetSize(); Int_t volId = volIds->At(0); @@ -161,7 +192,6 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, "fPoints.="<0){ + if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()MakeResiduals(); - AliTrackPoint p2local; for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){ - Int_t index = pindex[ipoint]; - fPoints->GetPoint(p,index); - if (GetPCA(p,p2)) { - plocal = p.Rotate(fAlpha); - // plocal.Rotate(fAlpha); - Float_t xyz[3],xyz2[3]; - p2local = plocal; - plocal.GetXYZ(xyz); - xyz2[0] = xyz[0]; - xyz2[1] = GetYat(xyz[0]); - xyz2[2] = GetZat(xyz[0]); - p2local.SetXYZ(xyz2); - lPVolId->AddPoint(ipoint,&plocal); - lPTrack->AddPoint(ipoint,&p2local); - } + 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[3]; + 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 @@ -290,17 +325,23 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, Int_t volId = volIds->At(0); Int_t modId =0; Int_t layer = AliAlignObj::VolUIDToLayer(volId,modId); + Int_t volIdFit = volIdsFit->At(0); + Int_t modIdFit =0; + Int_t layerFit = AliAlignObj::VolUIDToLayer(volIdFit,modIdFit); (*fDebugStream)<<"Fit"<< "VolId="<GetN(); fConv = kTRUE; } + // + // + // + 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); } + + //_____________________________________________________________________________ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const { @@ -400,6 +461,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+ @@ -416,13 +483,66 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c newcov[5] = cov[5]; p2.SetXYZ(xsecond,ysecond,zsecond,newcov); - if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){ + 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); +} + diff --git a/STEER/AliTrackFitterRieman.h b/STEER/AliTrackFitterRieman.h index 3ec9ffdaff4..9f4ac96ed65 100644 --- a/STEER/AliTrackFitterRieman.h +++ b/STEER/AliTrackFitterRieman.h @@ -33,25 +33,31 @@ class AliTrackFitterRieman : public AliTrackFitter{ Bool_t GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const; void SetMaxDelta(Float_t maxDelta) { fMaxDelta = maxDelta;} Float_t GetMaxDelta() const { return fMaxDelta;} - + void SetCorrection(Bool_t correction){ fBCorrection=correction;} + Bool_t GetCorrection() const {return fBCorrection ;} void Reset(); void AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz); void Update(); Double_t GetC() const {return fRieman->GetC();} - Double_t GetYat(Double_t x) const {return fRieman->GetYat(x);} - Double_t GetZat(Double_t x) const {return fRieman->GetZat(x);} + Double_t GetYat(Double_t x) const; + Double_t GetZat(Double_t x) const; Double_t GetDYat(Double_t x) const {return fRieman->GetDYat(x);} - Double_t GetDZat(Double_t x) const {return fRieman->GetDZat(x);} + Double_t GetDZat(Double_t x) const {return fRieman->GetDZat(x);} + Double_t GetErrY2at(Double_t x) const; + Double_t GetErrZ2at(Double_t x) const; + Bool_t GetXYZat(Double_t r, Float_t *xyz) const {return fRieman->GetXYZat(r, fAlpha,xyz);} AliRieman *GetRieman() const {return fRieman;} protected: - + Bool_t fBCorrection; //add correction for non-helicity Double_t fAlpha; //angle to transform to the fitting coordinate system Int_t fNUsed; //actual number of space-points used in the fit Bool_t fConv; //indicates convergation Float_t fMaxDelta; // maximal allowed delta in PCA exported for PCA minimization AliRieman *fRieman; // rieman fitter + Double_t fCorrY[4]; // correction polynom coef + Double_t fCorrZ[4]; // correction polynom coef private: TTreeSRedirector *fDebugStream; //!debug streamer ClassDef(AliTrackFitterRieman,2) // Fast fit of helices on ITS RecPoints diff --git a/STEER/AliTrackResidualsLinear.cxx b/STEER/AliTrackResidualsLinear.cxx index 840bcb6774d..7f34d031188 100644 --- a/STEER/AliTrackResidualsLinear.cxx +++ b/STEER/AliTrackResidualsLinear.cxx @@ -53,6 +53,7 @@ AliTrackResidualsLinear::AliTrackResidualsLinear(): fFixed[ipar] = 0;; fParams[ipar] = 0; } + for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;} } //______________________________________________________________________________ @@ -68,6 +69,7 @@ AliTrackResidualsLinear::AliTrackResidualsLinear(Int_t ntracks): fFixed[ipar] = 0; fParams[ipar] = 0; } + for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;} } //______________________________________________________________________________ @@ -84,6 +86,8 @@ AliTrackResidualsLinear::AliTrackResidualsLinear(const AliTrackResidualsLinear & fFixed[ipar] = res.fFixed[ipar]; fParams[ipar] = res.fParams[ipar]; } + for (Int_t icov=0; icov<36; icov++){ fCovar[icov]= res.fCovar[icov];} + fChi2Orig = res.fChi2Orig; } //______________________________________________________________________________ @@ -154,7 +158,7 @@ void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime) mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4]; mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5]; mcov+=mcovp; - mcov.Invert(); + //mcov.Invert(); if (!mcov.IsValid()) return; TMatrixD mcovBack = mcov; // for debug purposes // @@ -184,7 +188,7 @@ void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime) TMatrixD deltaT(matrixV, TMatrixD::kTransposeMult, deltaR); // tranformed delta TMatrixD mparamT(matrixV,TMatrixD::kTransposeMult, mparam); // tranformed linear transformation - if (0){ + if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>2){ // // debug part // @@ -210,20 +214,50 @@ void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime) printf("Rotated param matrix\n"); mparamT.Print(); // + // + printf("Trans Matrix:\n"); + matrixV.Print(); + printf("Delta Orig\n"); + deltaR.Print(); + printf("Delta Rotated"); + deltaT.Print(); + // + // } // - for (Int_t idim = 0; idim<3; idim++){ + Double_t sumChi2=0; + for (Int_t idim = 1; idim<3; idim++){ Double_t yf; // input points to fit in TLinear fitter Double_t xf[6]; // input points to fit yf = deltaT(idim,0); for (Int_t ipar =0; ipar<6; ipar++) xf[ipar] = mparamT(idim,ipar); if (covDiagonal[idim]>0.){ - fFitter->AddPoint(xf,yf, TMath::Sqrt(1/covDiagonal[idim])); + fFitter->AddPoint(xf,yf, TMath::Sqrt(covDiagonal[idim])); + // accumulate chi2 + Double_t chi2 = (yf*yf)/covDiagonal[idim]; + fChi2Orig += (yf*yf)/covDiagonal[idim]; + if (chi2>100 && AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){ + printf("Too big chi2- %f\n",chi2); + printf("Delta Orig\n"); + deltaR.Print(); + printf("Delta Rotated"); + deltaT.Print(); + matrixV.Print(); + printf("Too big chi2 - End\n"); + } + sumChi2+=chi2; + } + else{ + printf("Bug\n"); } - // accumulate chi2 - fChi2Orig += (yf*yf)*covDiagonal[idim]; } - fNdf +=3; + if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){ + TMatrixD matChi0=(mcov.Invert()*deltaR); + TMatrixD matChi2=deltaR.T()*matChi0; + printf("Chi2:\t%f\t%f", matChi2(0,0), sumChi2); + } + + fNdf +=2; } //______________________________________________________________________________ @@ -248,7 +282,7 @@ Bool_t AliTrackResidualsLinear::Update() } // fFitter->ReleaseParameter(0); - for (Int_t ipar=0; ipar<7; ipar++) { + for (Int_t ipar=0; ipar<6; ipar++) { if (fBFixed[ipar]) fFitter->ReleaseParameter(ipar+1); } @@ -264,6 +298,12 @@ Bool_t AliTrackResidualsLinear::Update() fParams[3] = vector[4]; fParams[4] = vector[5]; fParams[5] = vector[6]; + TMatrixD covar(7,7); + fFitter->GetCovarianceMatrix(covar); + for (Int_t i0=0; i0 <6; i0++) + for (Int_t j0=0; j0 <6; j0++){ + fCovar[i0*6+j0] = covar(i0+1,j0+1); + } // fAlignObj->SetPars(fParams[0], fParams[1], fParams[2], TMath::RadToDeg()*fParams[3], diff --git a/STEER/AliTrackResidualsLinear.h b/STEER/AliTrackResidualsLinear.h index 227c2497cea..18b3a56746b 100644 --- a/STEER/AliTrackResidualsLinear.h +++ b/STEER/AliTrackResidualsLinear.h @@ -25,13 +25,16 @@ class AliTrackResidualsLinear : public AliTrackResiduals { Bool_t Minimize(); void SetRobust(Float_t fraction){fFraction=fraction;} void FixParameter(Int_t par, Float_t value) {fBFixed[par]=kTRUE; fFixed[par]= value;} + const Double_t * GetParameters() const { return fParams;} + const Double_t * GetCovariance() const { return fCovar;} void ReleaseParameter(Int_t par) {fBFixed[par]=kFALSE;} protected: Bool_t Update(); void AddPoints(AliTrackPoint &p, AliTrackPoint &pprime); TLinearFitter *fFitter; // !Linear fitter Float_t fFraction; // fraction of points for robust fitter if less 0 - no robust fitter invoked - Float_t fParams[6]; // resulting parameters + Double_t fParams[6]; // resulting parameters + Double_t fCovar[36]; // covariance matrix Float_t fFixed[6]; // the fixed values of parameters Float_t fBFixed[6]; // the flag for fixing parameter Double_t fChi2Orig; // original chi2 -- 2.43.0