]>
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++){ |
bf9a6ef9 | 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++){ |
bf9a6ef9 | 68 | fParams[ipar] = 0; |
69 | } | |
dfce6628 | 70 | for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;} |
bf9a6ef9 | 71 | } |
72 | ||
73 | //______________________________________________________________________________ | |
74 | AliTrackResidualsLinear::AliTrackResidualsLinear(const AliTrackResidualsLinear &res): | |
75e3794b | 75 | AliTrackResiduals(res), |
76 | fFitter(new TLinearFitter(*(res.fFitter))), | |
77 | fFraction(res.fFraction), | |
78 | fChi2Orig(res.fChi2Orig) | |
bf9a6ef9 | 79 | { |
80 | // Copy constructor | |
81 | //.. | |
bf9a6ef9 | 82 | for (Int_t ipar=0; ipar<6; ipar++){ |
bf9a6ef9 | 83 | fParams[ipar] = res.fParams[ipar]; |
84 | } | |
dfce6628 | 85 | for (Int_t icov=0; icov<36; icov++){ fCovar[icov]= res.fCovar[icov];} |
86 | fChi2Orig = res.fChi2Orig; | |
bf9a6ef9 | 87 | } |
88 | ||
89 | //______________________________________________________________________________ | |
90 | AliTrackResidualsLinear &AliTrackResidualsLinear::operator= (const AliTrackResidualsLinear& res) | |
91 | { | |
92 | // Assignment operator | |
93 | ((AliTrackResiduals *)this)->operator=(res); | |
94 | return *this; | |
95 | } | |
96 | //______________________________________________________________________________ | |
97 | AliTrackResidualsLinear::~AliTrackResidualsLinear() | |
98 | { | |
99 | // | |
100 | // | |
101 | // | |
102 | delete fFitter; | |
103 | } | |
104 | ||
105 | ||
106 | //______________________________________________________________________________ | |
107 | Bool_t AliTrackResidualsLinear::Minimize() | |
108 | { | |
109 | // Implementation of fast linear Chi2 minimizer | |
110 | // based on TLinear fitter | |
111 | // | |
112 | if (!fFitter) fFitter = new TLinearFitter(6,"hyp6"); | |
113 | fFitter->StoreData(kTRUE); | |
114 | fFitter->ClearPoints(); | |
115 | fChi2Orig = 0; | |
116 | AliTrackPoint p1,p2; | |
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); | |
122 | AddPoints(p1,p2); | |
123 | } | |
124 | } | |
125 | Bool_t isOK = Update(); | |
126 | if (!isOK) return isOK; | |
127 | // | |
128 | TGeoHMatrix matrix; | |
129 | fAlignObj->GetMatrix(matrix); | |
130 | return isOK; | |
131 | } | |
132 | ||
133 | ||
134 | ||
135 | //______________________________________________________________________________ | |
136 | void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime) | |
137 | { | |
138 | // | |
139 | // add points to linear fitter - option with correlation betwee measurement in different dimensions | |
140 | // p1 - local point | |
141 | // pprime - track extrapolation point | |
142 | // | |
143 | Float_t xyz[3],xyzp[3]; | |
144 | Float_t cov[6],covp[6]; | |
145 | p.GetXYZ(xyz,cov); pprime.GetXYZ(xyzp,covp); | |
146 | // | |
147 | // | |
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]; | |
156 | mcov+=mcovp; | |
dfce6628 | 157 | //mcov.Invert(); |
bf9a6ef9 | 158 | if (!mcov.IsValid()) return; |
159 | TMatrixD mcovBack = mcov; // for debug purposes | |
160 | // | |
161 | // decompose matrix | |
162 | // | |
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 | |
167 | // | |
168 | // residual vector | |
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]); | |
173 | // | |
174 | // parametrization matrix | |
175 | // | |
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 | |
183 | // | |
184 | ||
185 | TMatrixD deltaT(matrixV, TMatrixD::kTransposeMult, deltaR); // tranformed delta | |
186 | TMatrixD mparamT(matrixV,TMatrixD::kTransposeMult, mparam); // tranformed linear transformation | |
dfce6628 | 187 | if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>2){ |
bf9a6ef9 | 188 | // |
189 | // debug part | |
190 | // | |
191 | // covDiag = U^-1 * mcov * V -- diagonalization of covariance matrix | |
192 | ||
193 | TMatrixD matrixU = svd.GetU(); // transformation matrix to diagonalize covariance matrix | |
194 | TMatrixD matrixUI= svd.GetU(); | |
195 | matrixUI.Invert(); | |
196 | // | |
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"); | |
201 | test2.Print(); | |
202 | printf("Test matrix 0 - should be unit\n"); | |
203 | test0.Print(); | |
204 | printf("Test matrix 1 - should be diagonal\n"); | |
205 | test1.Print(); | |
206 | printf("Diagonal matrix\n"); | |
207 | svd.GetSig().Print(); | |
208 | printf("Original param matrix\n"); | |
209 | mparam.Print(); | |
210 | printf("Rotated param matrix\n"); | |
211 | mparamT.Print(); | |
212 | // | |
dfce6628 | 213 | // |
214 | printf("Trans Matrix:\n"); | |
215 | matrixV.Print(); | |
216 | printf("Delta Orig\n"); | |
217 | deltaR.Print(); | |
218 | printf("Delta Rotated"); | |
219 | deltaT.Print(); | |
220 | // | |
221 | // | |
bf9a6ef9 | 222 | } |
223 | // | |
dfce6628 | 224 | Double_t sumChi2=0; |
225 | for (Int_t idim = 1; idim<3; idim++){ | |
bf9a6ef9 | 226 | Double_t yf; // input points to fit in TLinear fitter |
227 | Double_t xf[6]; // input points to fit | |
228 | yf = deltaT(idim,0); | |
229 | for (Int_t ipar =0; ipar<6; ipar++) xf[ipar] = mparamT(idim,ipar); | |
230 | if (covDiagonal[idim]>0.){ | |
dfce6628 | 231 | fFitter->AddPoint(xf,yf, TMath::Sqrt(covDiagonal[idim])); |
232 | // accumulate chi2 | |
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"); | |
238 | deltaR.Print(); | |
239 | printf("Delta Rotated"); | |
240 | deltaT.Print(); | |
241 | matrixV.Print(); | |
242 | printf("Too big chi2 - End\n"); | |
243 | } | |
244 | sumChi2+=chi2; | |
245 | } | |
246 | else{ | |
247 | printf("Bug\n"); | |
bf9a6ef9 | 248 | } |
bf9a6ef9 | 249 | } |
dfce6628 | 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); | |
254 | } | |
255 | ||
256 | fNdf +=2; | |
bf9a6ef9 | 257 | } |
258 | ||
259 | //______________________________________________________________________________ | |
260 | Bool_t AliTrackResidualsLinear::Update() | |
261 | { | |
262 | // Find the alignment parameters | |
263 | // using TLinear fitter + fill data containers | |
264 | // | |
265 | // | |
266 | fFitter->Eval(); | |
267 | // | |
268 | // TLinear fitter put as first parameter offset - fixing parameter shifted by one | |
269 | // | |
270 | fFitter->FixParameter(0); | |
271 | for (Int_t ipar =0; ipar<6; ipar++){ | |
272 | if (fBFixed[ipar]) fFitter->FixParameter(ipar+1,fFixed[ipar]); | |
273 | } | |
274 | if (fFraction>0.5) { | |
275 | fFitter->EvalRobust(fFraction); | |
276 | }else{ | |
277 | fFitter->Eval(); | |
278 | } | |
279 | // | |
280 | fFitter->ReleaseParameter(0); | |
dfce6628 | 281 | for (Int_t ipar=0; ipar<6; ipar++) { |
bf9a6ef9 | 282 | if (fBFixed[ipar]) fFitter->ReleaseParameter(ipar+1); |
283 | } | |
284 | ||
285 | ||
286 | // | |
287 | fChi2 = fFitter->GetChisquare(); | |
288 | fNdf -= 6; | |
289 | TVectorD vector(7); | |
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]; | |
dfce6628 | 297 | TMatrixD covar(7,7); |
298 | fFitter->GetCovarianceMatrix(covar); | |
299 | for (Int_t i0=0; i0 <6; i0++) | |
300 | for (Int_t j0=0; j0 <6; j0++){ | |
301 | fCovar[i0*6+j0] = covar(i0+1,j0+1); | |
302 | } | |
bf9a6ef9 | 303 | // |
304 | fAlignObj->SetPars(fParams[0], fParams[1], fParams[2], | |
305 | TMath::RadToDeg()*fParams[3], | |
306 | TMath::RadToDeg()*fParams[4], | |
307 | TMath::RadToDeg()*fParams[5]); | |
308 | return kTRUE; | |
309 | } |