1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-----------------------------------------------------------------
19 // Implementation of the derived class for track residuals
20 // based on linear chi2 minimization
21 // The minimization relies on the fact that the alignment parameters
22 // (angles and translations) are small.
23 // TLinearFitter used for minimization
24 // Possibility to fix Paramaters
25 // FixParameter()ReleaseParameter();
26 // Possibility to define fraction of outliers to be skipped
28 // marian.ivanov@cern.ch
30 //-----------------------------------------------------------------
33 #include <TGeoMatrix.h>
36 #include "AliAlignObj.h"
37 #include "AliTrackPointArray.h"
38 #include "AliTrackResidualsLinear.h"
39 #include "AliAlignObj.h"
40 #include "TLinearFitter.h"
41 #include "TDecompSVD.h"
43 ClassImp(AliTrackResidualsLinear)
45 //______________________________________________________________________________
46 AliTrackResidualsLinear::AliTrackResidualsLinear():
52 // Default constructor
53 for (Int_t ipar=0; ipar<6; ipar++){
56 for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;}
59 //______________________________________________________________________________
60 AliTrackResidualsLinear::AliTrackResidualsLinear(Int_t ntracks):
61 AliTrackResiduals(ntracks),
62 fFitter(new TLinearFitter(6,"hyp6")),
67 for (Int_t ipar=0; ipar<6; ipar++){
70 for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;}
73 //______________________________________________________________________________
74 AliTrackResidualsLinear::AliTrackResidualsLinear(const AliTrackResidualsLinear &res):
75 AliTrackResiduals(res),
76 fFitter(new TLinearFitter(*(res.fFitter))),
77 fFraction(res.fFraction),
78 fChi2Orig(res.fChi2Orig)
82 for (Int_t ipar=0; ipar<6; ipar++){
83 fParams[ipar] = res.fParams[ipar];
85 for (Int_t icov=0; icov<36; icov++){ fCovar[icov]= res.fCovar[icov];}
86 fChi2Orig = res.fChi2Orig;
89 //______________________________________________________________________________
90 AliTrackResidualsLinear &AliTrackResidualsLinear::operator= (const AliTrackResidualsLinear& res)
92 // Assignment operator
93 ((AliTrackResiduals *)this)->operator=(res);
96 //______________________________________________________________________________
97 AliTrackResidualsLinear::~AliTrackResidualsLinear()
106 //______________________________________________________________________________
107 Bool_t AliTrackResidualsLinear::Minimize()
109 // Implementation of fast linear Chi2 minimizer
110 // based on TLinear fitter
112 if (!fFitter) fFitter = new TLinearFitter(6,"hyp6");
113 fFitter->StoreData(kTRUE);
114 fFitter->ClearPoints();
117 for (Int_t itrack = 0; itrack < fLast; itrack++) {
118 if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
119 for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
120 fVolArray[itrack]->GetPoint(p1,ipoint);
121 fTrackArray[itrack]->GetPoint(p2,ipoint);
125 Bool_t isOK = Update();
126 if (!isOK) return isOK;
129 fAlignObj->GetMatrix(matrix);
135 //______________________________________________________________________________
136 void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime)
139 // add points to linear fitter - option with correlation betwee measurement in different dimensions
141 // pprime - track extrapolation point
143 Float_t xyz[3],xyzp[3];
144 Float_t cov[6],covp[6];
145 p.GetXYZ(xyz,cov); pprime.GetXYZ(xyzp,covp);
148 TMatrixD mcov(3,3); // local point covariance
149 mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
150 mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
151 mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
152 TMatrixD mcovp(3,3); // extrapolation point covariance
153 mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
154 mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
155 mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
158 if (!mcov.IsValid()) return;
159 TMatrixD mcovBack = mcov; // for debug purposes
163 TDecompSVD svd(mcov); // mcov = svd.fU * covDiagonal * svd.fV.Invert
164 if (!svd.Decompose()) return; // decomposition failed
165 TMatrixD matrixV = svd.GetV(); // transformation matrix to diagonalize covariance matrix
166 Double_t covDiagonal[3] = {svd.GetSig()[0],svd.GetSig()[1],svd.GetSig()[2]}; // diagonalized covariance matrix
169 TMatrixD deltaR(3,1);
170 deltaR(0,0) = (xyzp[0]-xyz[0]);
171 deltaR(1,0) = (xyzp[1]-xyz[1]);
172 deltaR(2,0) = (xyzp[2]-xyz[2]);
174 // parametrization matrix
176 TMatrixD mparam(3,6);
177 mparam(0,0) = 1; mparam(1,0) = 0; mparam(2,0) = 0; // xshift
178 mparam(0,1) = 0; mparam(1,1) = 1; mparam(2,1) = 0; // yshift
179 mparam(0,2) = 0; mparam(1,2) = 0; mparam(2,2) = 1; // zshift
180 mparam(0,3) = 0; mparam(1,3) =-xyz[2]; mparam(2,3) = xyz[1]; // x rotation
181 mparam(0,4) = xyz[2]; mparam(1,4) = 0; mparam(2,4) =-xyz[0]; // y rotation
182 mparam(0,5) =-xyz[1]; mparam(1,5) = xyz[0]; mparam(2,5) = 0; // z rotation
185 TMatrixD deltaT(matrixV, TMatrixD::kTransposeMult, deltaR); // tranformed delta
186 TMatrixD mparamT(matrixV,TMatrixD::kTransposeMult, mparam); // tranformed linear transformation
187 if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>2){
191 // covDiag = U^-1 * mcov * V -- diagonalization of covariance matrix
193 TMatrixD matrixU = svd.GetU(); // transformation matrix to diagonalize covariance matrix
194 TMatrixD matrixUI= svd.GetU();
197 TMatrixD test0 = matrixUI*matrixV; // test matrix - should be unit matrix
198 TMatrixD test1 = matrixUI*mcovBack*matrixV; // test matrix - diagonal - should be diagonal with covDiagonal on diag
199 TMatrixD test2 = matrixU.T()*matrixV; // test ortogonality - shoul be unit
200 printf("Test matrix 2 - should be unit\n");
202 printf("Test matrix 0 - should be unit\n");
204 printf("Test matrix 1 - should be diagonal\n");
206 printf("Diagonal matrix\n");
207 svd.GetSig().Print();
208 printf("Original param matrix\n");
210 printf("Rotated param matrix\n");
214 printf("Trans Matrix:\n");
216 printf("Delta Orig\n");
218 printf("Delta Rotated");
225 for (Int_t idim = 1; idim<3; idim++){
226 Double_t yf; // input points to fit in TLinear fitter
227 Double_t xf[6]; // input points to fit
229 for (Int_t ipar =0; ipar<6; ipar++) xf[ipar] = mparamT(idim,ipar);
230 if (covDiagonal[idim]>0.){
231 fFitter->AddPoint(xf,yf, TMath::Sqrt(covDiagonal[idim]));
233 Double_t chi2 = (yf*yf)/covDiagonal[idim];
234 fChi2Orig += (yf*yf)/covDiagonal[idim];
235 if (chi2>100 && AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){
236 printf("Too big chi2- %f\n",chi2);
237 printf("Delta Orig\n");
239 printf("Delta Rotated");
242 printf("Too big chi2 - End\n");
250 if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){
251 TMatrixD matChi0=(mcov.Invert()*deltaR);
252 TMatrixD matChi2=deltaR.T()*matChi0;
253 printf("Chi2:\t%f\t%f", matChi2(0,0), sumChi2);
259 //______________________________________________________________________________
260 Bool_t AliTrackResidualsLinear::Update()
262 // Find the alignment parameters
263 // using TLinear fitter + fill data containers
268 // TLinear fitter put as first parameter offset - fixing parameter shifted by one
270 fFitter->FixParameter(0);
271 for (Int_t ipar =0; ipar<6; ipar++){
272 if (fBFixed[ipar]) fFitter->FixParameter(ipar+1,fFixed[ipar]);
275 fFitter->EvalRobust(fFraction);
280 fFitter->ReleaseParameter(0);
281 for (Int_t ipar=0; ipar<6; ipar++) {
282 if (fBFixed[ipar]) fFitter->ReleaseParameter(ipar+1);
287 fChi2 = fFitter->GetChisquare();
290 fFitter->GetParameters(vector);
291 fParams[0] = vector[1];
292 fParams[1] = vector[2];
293 fParams[2] = vector[3];
294 fParams[3] = vector[4];
295 fParams[4] = vector[5];
296 fParams[5] = vector[6];
298 fFitter->GetCovarianceMatrix(covar);
300 for (Int_t i0=0; i0 <6; i0++){
301 for (Int_t j0=0; j0 <6; j0++){
302 fCovar[i0*6+j0] = covar(i0+1,j0+1);
305 Double_t covmatrarray[21];
306 for(Int_t j=0;j<6;j++){
307 for(Int_t i=j;i<6;i++){
308 covmatrarray[i*(i+1)/2+j]=fCovar[i+6*j];
313 fAlignObj->SetPars(fParams[0], fParams[1], fParams[2],
314 TMath::RadToDeg()*fParams[3],
315 TMath::RadToDeg()*fParams[4],
316 TMath::RadToDeg()*fParams[5]);
317 fAlignObj->SetCorrMatrix(covmatrarray);