]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliTrackResidualsLinear.cxx
moved AliSelector, AliSelectorRL to STEER
[u/mrichter/AliRoot.git] / STEER / AliTrackResidualsLinear.cxx
CommitLineData
bf9a6ef9 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16//-----------------------------------------------------------------
17// Implementation of the derived class for track residuals
18// based on linear chi2 minimization
19// The minimization relies on the fact that the alignment parameters
20// (angles and translations) are small.
21// TLinearFitter used for minimization
22// Possibility to fix Paramaters
23// FixParameter()ReleaseParameter();
24// Possibility to define fraction of outliers to be skipped
25//
26// marian.ivanov@cern.ch
27//
28//-----------------------------------------------------------------
29
30#include <TMinuit.h>
31#include <TGeoMatrix.h>
32
33#include "AliLog.h"
34#include "AliAlignObj.h"
35#include "AliTrackPointArray.h"
36#include "AliTrackResidualsLinear.h"
37#include "AliAlignObj.h"
38#include "TLinearFitter.h"
39#include "TDecompSVD.h"
40
41ClassImp(AliTrackResidualsLinear)
42
43//______________________________________________________________________________
75e3794b 44AliTrackResidualsLinear::AliTrackResidualsLinear():
45 AliTrackResiduals(),
46 fFitter(0),
47 fFraction(-1),
48 fChi2Orig(0)
bf9a6ef9 49{
50 // Default constructor
bf9a6ef9 51 for (Int_t ipar=0; ipar<6; ipar++){
52 fBFixed[ipar] = kFALSE;
53 fFixed[ipar] = 0;;
54 fParams[ipar] = 0;
55 }
dfce6628 56 for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;}
bf9a6ef9 57}
58
59//______________________________________________________________________________
60AliTrackResidualsLinear::AliTrackResidualsLinear(Int_t ntracks):
75e3794b 61 AliTrackResiduals(ntracks),
62 fFitter(new TLinearFitter(6,"hyp6")),
63 fFraction(-1),
64 fChi2Orig(0)
bf9a6ef9 65{
66 // Constructor
bf9a6ef9 67 for (Int_t ipar=0; ipar<6; ipar++){
68 fBFixed[ipar] = kFALSE;
69 fFixed[ipar] = 0;
70 fParams[ipar] = 0;
71 }
dfce6628 72 for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;}
bf9a6ef9 73}
74
75//______________________________________________________________________________
76AliTrackResidualsLinear::AliTrackResidualsLinear(const AliTrackResidualsLinear &res):
75e3794b 77 AliTrackResiduals(res),
78 fFitter(new TLinearFitter(*(res.fFitter))),
79 fFraction(res.fFraction),
80 fChi2Orig(res.fChi2Orig)
bf9a6ef9 81{
82 // Copy constructor
83 //..
bf9a6ef9 84 for (Int_t ipar=0; ipar<6; ipar++){
85 fBFixed[ipar] = res.fBFixed[ipar];
86 fFixed[ipar] = res.fFixed[ipar];
87 fParams[ipar] = res.fParams[ipar];
88 }
dfce6628 89 for (Int_t icov=0; icov<36; icov++){ fCovar[icov]= res.fCovar[icov];}
90 fChi2Orig = res.fChi2Orig;
bf9a6ef9 91}
92
93//______________________________________________________________________________
94AliTrackResidualsLinear &AliTrackResidualsLinear::operator= (const AliTrackResidualsLinear& res)
95{
96 // Assignment operator
97 ((AliTrackResiduals *)this)->operator=(res);
98 return *this;
99}
100//______________________________________________________________________________
101AliTrackResidualsLinear::~AliTrackResidualsLinear()
102{
103 //
104 //
105 //
106 delete fFitter;
107}
108
109
110//______________________________________________________________________________
111Bool_t AliTrackResidualsLinear::Minimize()
112{
113 // Implementation of fast linear Chi2 minimizer
114 // based on TLinear fitter
115 //
116 if (!fFitter) fFitter = new TLinearFitter(6,"hyp6");
117 fFitter->StoreData(kTRUE);
118 fFitter->ClearPoints();
119 fChi2Orig = 0;
120 AliTrackPoint p1,p2;
121 for (Int_t itrack = 0; itrack < fLast; itrack++) {
122 if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
123 for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
124 fVolArray[itrack]->GetPoint(p1,ipoint);
125 fTrackArray[itrack]->GetPoint(p2,ipoint);
126 AddPoints(p1,p2);
127 }
128 }
129 Bool_t isOK = Update();
130 if (!isOK) return isOK;
131 //
132 TGeoHMatrix matrix;
133 fAlignObj->GetMatrix(matrix);
134 return isOK;
135}
136
137
138
139//______________________________________________________________________________
140void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime)
141{
142 //
143 // add points to linear fitter - option with correlation betwee measurement in different dimensions
144 // p1 - local point
145 // pprime - track extrapolation point
146 //
147 Float_t xyz[3],xyzp[3];
148 Float_t cov[6],covp[6];
149 p.GetXYZ(xyz,cov); pprime.GetXYZ(xyzp,covp);
150 //
151 //
152 TMatrixD mcov(3,3); // local point covariance
153 mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
154 mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
155 mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
156 TMatrixD mcovp(3,3); // extrapolation point covariance
157 mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
158 mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
159 mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
160 mcov+=mcovp;
dfce6628 161 //mcov.Invert();
bf9a6ef9 162 if (!mcov.IsValid()) return;
163 TMatrixD mcovBack = mcov; // for debug purposes
164 //
165 // decompose matrix
166 //
167 TDecompSVD svd(mcov); // mcov = svd.fU * covDiagonal * svd.fV.Invert
168 if (!svd.Decompose()) return; // decomposition failed
169 TMatrixD matrixV = svd.GetV(); // transformation matrix to diagonalize covariance matrix
170 Double_t covDiagonal[3] = {svd.GetSig()[0],svd.GetSig()[1],svd.GetSig()[2]}; // diagonalized covariance matrix
171 //
172 // residual vector
173 TMatrixD deltaR(3,1);
174 deltaR(0,0) = (xyzp[0]-xyz[0]);
175 deltaR(1,0) = (xyzp[1]-xyz[1]);
176 deltaR(2,0) = (xyzp[2]-xyz[2]);
177 //
178 // parametrization matrix
179 //
180 TMatrixD mparam(3,6);
181 mparam(0,0) = 1; mparam(1,0) = 0; mparam(2,0) = 0; // xshift
182 mparam(0,1) = 0; mparam(1,1) = 1; mparam(2,1) = 0; // yshift
183 mparam(0,2) = 0; mparam(1,2) = 0; mparam(2,2) = 1; // zshift
184 mparam(0,3) = 0; mparam(1,3) =-xyz[2]; mparam(2,3) = xyz[1]; // x rotation
185 mparam(0,4) = xyz[2]; mparam(1,4) = 0; mparam(2,4) =-xyz[0]; // y rotation
186 mparam(0,5) =-xyz[1]; mparam(1,5) = xyz[0]; mparam(2,5) = 0; // z rotation
187 //
188
189 TMatrixD deltaT(matrixV, TMatrixD::kTransposeMult, deltaR); // tranformed delta
190 TMatrixD mparamT(matrixV,TMatrixD::kTransposeMult, mparam); // tranformed linear transformation
dfce6628 191 if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>2){
bf9a6ef9 192 //
193 // debug part
194 //
195 // covDiag = U^-1 * mcov * V -- diagonalization of covariance matrix
196
197 TMatrixD matrixU = svd.GetU(); // transformation matrix to diagonalize covariance matrix
198 TMatrixD matrixUI= svd.GetU();
199 matrixUI.Invert();
200 //
201 TMatrixD test0 = matrixUI*matrixV; // test matrix - should be unit matrix
202 TMatrixD test1 = matrixUI*mcovBack*matrixV; // test matrix - diagonal - should be diagonal with covDiagonal on diag
203 TMatrixD test2 = matrixU.T()*matrixV; // test ortogonality - shoul be unit
204 printf("Test matrix 2 - should be unit\n");
205 test2.Print();
206 printf("Test matrix 0 - should be unit\n");
207 test0.Print();
208 printf("Test matrix 1 - should be diagonal\n");
209 test1.Print();
210 printf("Diagonal matrix\n");
211 svd.GetSig().Print();
212 printf("Original param matrix\n");
213 mparam.Print();
214 printf("Rotated param matrix\n");
215 mparamT.Print();
216 //
dfce6628 217 //
218 printf("Trans Matrix:\n");
219 matrixV.Print();
220 printf("Delta Orig\n");
221 deltaR.Print();
222 printf("Delta Rotated");
223 deltaT.Print();
224 //
225 //
bf9a6ef9 226 }
227 //
dfce6628 228 Double_t sumChi2=0;
229 for (Int_t idim = 1; idim<3; idim++){
bf9a6ef9 230 Double_t yf; // input points to fit in TLinear fitter
231 Double_t xf[6]; // input points to fit
232 yf = deltaT(idim,0);
233 for (Int_t ipar =0; ipar<6; ipar++) xf[ipar] = mparamT(idim,ipar);
234 if (covDiagonal[idim]>0.){
dfce6628 235 fFitter->AddPoint(xf,yf, TMath::Sqrt(covDiagonal[idim]));
236 // accumulate chi2
237 Double_t chi2 = (yf*yf)/covDiagonal[idim];
238 fChi2Orig += (yf*yf)/covDiagonal[idim];
239 if (chi2>100 && AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){
240 printf("Too big chi2- %f\n",chi2);
241 printf("Delta Orig\n");
242 deltaR.Print();
243 printf("Delta Rotated");
244 deltaT.Print();
245 matrixV.Print();
246 printf("Too big chi2 - End\n");
247 }
248 sumChi2+=chi2;
249 }
250 else{
251 printf("Bug\n");
bf9a6ef9 252 }
bf9a6ef9 253 }
dfce6628 254 if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){
255 TMatrixD matChi0=(mcov.Invert()*deltaR);
256 TMatrixD matChi2=deltaR.T()*matChi0;
257 printf("Chi2:\t%f\t%f", matChi2(0,0), sumChi2);
258 }
259
260 fNdf +=2;
bf9a6ef9 261}
262
263//______________________________________________________________________________
264Bool_t AliTrackResidualsLinear::Update()
265{
266 // Find the alignment parameters
267 // using TLinear fitter + fill data containers
268 //
269 //
270 fFitter->Eval();
271 //
272 // TLinear fitter put as first parameter offset - fixing parameter shifted by one
273 //
274 fFitter->FixParameter(0);
275 for (Int_t ipar =0; ipar<6; ipar++){
276 if (fBFixed[ipar]) fFitter->FixParameter(ipar+1,fFixed[ipar]);
277 }
278 if (fFraction>0.5) {
279 fFitter->EvalRobust(fFraction);
280 }else{
281 fFitter->Eval();
282 }
283 //
284 fFitter->ReleaseParameter(0);
dfce6628 285 for (Int_t ipar=0; ipar<6; ipar++) {
bf9a6ef9 286 if (fBFixed[ipar]) fFitter->ReleaseParameter(ipar+1);
287 }
288
289
290 //
291 fChi2 = fFitter->GetChisquare();
292 fNdf -= 6;
293 TVectorD vector(7);
294 fFitter->GetParameters(vector);
295 fParams[0] = vector[1];
296 fParams[1] = vector[2];
297 fParams[2] = vector[3];
298 fParams[3] = vector[4];
299 fParams[4] = vector[5];
300 fParams[5] = vector[6];
dfce6628 301 TMatrixD covar(7,7);
302 fFitter->GetCovarianceMatrix(covar);
303 for (Int_t i0=0; i0 <6; i0++)
304 for (Int_t j0=0; j0 <6; j0++){
305 fCovar[i0*6+j0] = covar(i0+1,j0+1);
306 }
bf9a6ef9 307 //
308 fAlignObj->SetPars(fParams[0], fParams[1], fParams[2],
309 TMath::RadToDeg()*fParams[3],
310 TMath::RadToDeg()*fParams[4],
311 TMath::RadToDeg()*fParams[5]);
312 return kTRUE;
313}