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