]>
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 | ||
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 | ||
41 | ClassImp(AliTrackResidualsLinear) | |
42 | ||
43 | //______________________________________________________________________________ | |
75e3794b | 44 | AliTrackResidualsLinear::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 | //______________________________________________________________________________ | |
60 | AliTrackResidualsLinear::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 | //______________________________________________________________________________ | |
76 | AliTrackResidualsLinear::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 | //______________________________________________________________________________ | |
94 | AliTrackResidualsLinear &AliTrackResidualsLinear::operator= (const AliTrackResidualsLinear& res) | |
95 | { | |
96 | // Assignment operator | |
97 | ((AliTrackResiduals *)this)->operator=(res); | |
98 | return *this; | |
99 | } | |
100 | //______________________________________________________________________________ | |
101 | AliTrackResidualsLinear::~AliTrackResidualsLinear() | |
102 | { | |
103 | // | |
104 | // | |
105 | // | |
106 | delete fFitter; | |
107 | } | |
108 | ||
109 | ||
110 | //______________________________________________________________________________ | |
111 | Bool_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 | //______________________________________________________________________________ | |
140 | void 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 | //______________________________________________________________________________ | |
264 | Bool_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 | } |