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