]>
Commit | Line | Data |
---|---|---|
46ae650f | 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 | ||
46ae650f | 18 | //----------------------------------------------------------------- |
19 | // Implementation of the derived class for track residuals | |
20 | // based on linear chi2 minimization (in approximation of | |
21 | // small alignment angles and translations) | |
22 | // | |
23 | //----------------------------------------------------------------- | |
24 | ||
3010c308 | 25 | #include <TMath.h> |
46ae650f | 26 | #include <TGeoMatrix.h> |
27 | ||
cc345ce3 | 28 | #include "AliLog.h" |
46ae650f | 29 | #include "AliAlignObj.h" |
30 | #include "AliTrackPointArray.h" | |
31 | #include "AliTrackResidualsFast.h" | |
32 | ||
4c84f475 | 33 | #include <TMatrixDSym.h> |
34 | #include <TMatrixDSymEigen.h> | |
35 | ||
46ae650f | 36 | ClassImp(AliTrackResidualsFast) |
37 | ||
38 | //______________________________________________________________________________ | |
75e3794b | 39 | AliTrackResidualsFast::AliTrackResidualsFast(): |
40 | AliTrackResiduals(), | |
41 | fSumR(0) | |
46ae650f | 42 | { |
43 | // Default constructor | |
cc345ce3 | 44 | for (Int_t i = 0; i < 27; i++) fSum[i] = 0; |
46ae650f | 45 | } |
46 | ||
47 | //______________________________________________________________________________ | |
48 | AliTrackResidualsFast::AliTrackResidualsFast(Int_t ntracks): | |
75e3794b | 49 | AliTrackResiduals(ntracks), |
50 | fSumR(0) | |
46ae650f | 51 | { |
52 | // Constructor | |
cc345ce3 | 53 | for (Int_t i = 0; i < 27; i++) fSum[i] = 0; |
46ae650f | 54 | } |
55 | ||
56 | //______________________________________________________________________________ | |
57 | AliTrackResidualsFast::AliTrackResidualsFast(const AliTrackResidualsFast &res): | |
75e3794b | 58 | AliTrackResiduals(res), |
59 | fSumR(res.fSumR) | |
46ae650f | 60 | { |
61 | // Copy constructor | |
cc345ce3 | 62 | for (Int_t i = 0; i < 27; i++) fSum[i] = res.fSum[i]; |
46ae650f | 63 | } |
64 | ||
65 | //______________________________________________________________________________ | |
66 | AliTrackResidualsFast &AliTrackResidualsFast::operator= (const AliTrackResidualsFast& res) | |
67 | { | |
68 | // Assignment operator | |
69 | ((AliTrackResiduals *)this)->operator=(res); | |
cc345ce3 | 70 | for (Int_t i = 0; i < 27; i++) fSum[i] = res.fSum[i]; |
46ae650f | 71 | fSumR = res.fSumR; |
72 | ||
73 | return *this; | |
74 | } | |
75 | ||
76 | //______________________________________________________________________________ | |
77 | Bool_t AliTrackResidualsFast::Minimize() | |
78 | { | |
79 | // Implementation of fast linear Chi2 | |
80 | // based minimization of track residuals sum | |
cc101660 | 81 | |
82 | // if(fBFixed[0]||fBFixed[1]||fBFixed[2]||fBFixed[3]||fBFixed[4]||fBFixed[5]) | |
83 | // AliError("Cannot yet fix parameters in this minimizer"); | |
29317b68 | 84 | |
85 | ||
cc345ce3 | 86 | for (Int_t i = 0; i < 27; i++) fSum[i] = 0; |
46ae650f | 87 | fSumR = 0; |
88 | ||
89 | AliTrackPoint p1,p2; | |
90 | ||
91 | for (Int_t itrack = 0; itrack < fLast; itrack++) { | |
92 | if (!fVolArray[itrack] || !fTrackArray[itrack]) continue; | |
93 | for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) { | |
94 | fVolArray[itrack]->GetPoint(p1,ipoint); | |
95 | fTrackArray[itrack]->GetPoint(p2,ipoint); | |
96 | AddPoints(p1,p2); | |
97 | } | |
98 | } | |
99 | ||
100 | return Update(); | |
101 | ||
102 | // debug info | |
103 | // Float_t chi2 = 0; | |
104 | // for (Int_t itrack = 0; itrack < fLast; itrack++) { | |
105 | // if (!fVolArray[itrack] || !fTrackArray[itrack]) continue; | |
106 | // for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) { | |
107 | // fVolArray[itrack]->GetPoint(p1,ipoint); | |
108 | // fAlignObj->Transform(p1); | |
109 | // fTrackArray[itrack]->GetPoint(p2,ipoint); | |
110 | // Float_t residual = p2.GetResidual(p1,kFALSE); | |
111 | // chi2 += residual; | |
112 | // } | |
113 | // } | |
114 | // printf("Final chi2 = %f\n",chi2); | |
115 | } | |
116 | ||
117 | //______________________________________________________________________________ | |
118 | void AliTrackResidualsFast::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime) | |
119 | { | |
120 | // Update the sums used for | |
121 | // the linear chi2 minimization | |
122 | Float_t xyz[3],xyzp[3]; | |
cc345ce3 | 123 | Float_t cov[6],covp[6]; |
124 | p.GetXYZ(xyz,cov); pprime.GetXYZ(xyzp,covp); | |
125 | TMatrixDSym mcov(3); | |
126 | mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2]; | |
127 | mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4]; | |
128 | mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5]; | |
129 | TMatrixDSym mcovp(3); | |
130 | mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2]; | |
131 | mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4]; | |
132 | mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5]; | |
133 | TMatrixDSym msum = mcov + mcovp; | |
4c84f475 | 134 | |
135 | ||
cc345ce3 | 136 | msum.Invert(); |
4c84f475 | 137 | |
138 | ||
cc345ce3 | 139 | if (!msum.IsValid()) return; |
140 | ||
141 | TMatrixD sums(3,1); | |
142 | sums(0,0) = (xyzp[0]-xyz[0]); | |
143 | sums(1,0) = (xyzp[1]-xyz[1]); | |
144 | sums(2,0) = (xyzp[2]-xyz[2]); | |
145 | TMatrixD sumst = sums.T(); sums.T(); | |
146 | ||
147 | TMatrixD mf(3,6); | |
148 | mf(0,0) = 1; mf(1,0) = 0; mf(2,0) = 0; | |
149 | mf(0,1) = 0; mf(1,1) = 1; mf(2,1) = 0; | |
150 | mf(0,2) = 0; mf(1,2) = 0; mf(2,2) = 1; | |
4c84f475 | 151 | mf(0,3) = 0; mf(1,3) = -xyz[2]; mf(2,3) = xyz[1]; |
cc345ce3 | 152 | mf(0,4) = xyz[2]; mf(1,4) = 0; mf(2,4) =-xyz[0]; |
153 | mf(0,5) =-xyz[1]; mf(1,5) = xyz[0]; mf(2,5) = 0; | |
cc101660 | 154 | |
155 | for(Int_t j=0;j<6;j++){ | |
156 | if(fBFixed[j]==kTRUE){ | |
157 | mf(0,j)=0.;mf(1,j)=0.;mf(2,j)=0.; | |
158 | } | |
159 | } | |
160 | ||
cc345ce3 | 161 | TMatrixD mft = mf.T(); mf.T(); |
162 | TMatrixD sums2 = mft * msum * sums; | |
163 | ||
164 | TMatrixD smatrix = mft * msum * mf; | |
165 | ||
166 | fSum[0] += smatrix(0,0); | |
167 | fSum[1] += smatrix(0,1); | |
168 | fSum[2] += smatrix(0,2); | |
169 | fSum[3] += smatrix(0,3); | |
170 | fSum[4] += smatrix(0,4); | |
171 | fSum[5] += smatrix(0,5); | |
172 | fSum[6] += smatrix(1,1); | |
173 | fSum[7] += smatrix(1,2); | |
174 | fSum[8] += smatrix(1,3); | |
175 | fSum[9] += smatrix(1,4); | |
176 | fSum[10]+= smatrix(1,5); | |
177 | fSum[11]+= smatrix(2,2); | |
178 | fSum[12]+= smatrix(2,3); | |
179 | fSum[13]+= smatrix(2,4); | |
180 | fSum[14]+= smatrix(2,5); | |
181 | fSum[15]+= smatrix(3,3); | |
182 | fSum[16]+= smatrix(3,4); | |
183 | fSum[17]+= smatrix(3,5); | |
184 | fSum[18]+= smatrix(4,4); | |
185 | fSum[19]+= smatrix(4,5); | |
186 | fSum[20]+= smatrix(5,5); | |
187 | fSum[21] += sums2(0,0); | |
188 | fSum[22] += sums2(1,0); | |
189 | fSum[23] += sums2(2,0); | |
190 | fSum[24] += sums2(3,0); | |
191 | fSum[25] += sums2(4,0); | |
192 | fSum[26] += sums2(5,0); | |
193 | ||
194 | TMatrixD tmp = sumst * msum * sums; | |
195 | fSumR += tmp(0,0); | |
46ae650f | 196 | |
197 | fNdf += 3; | |
198 | } | |
199 | ||
200 | //______________________________________________________________________________ | |
201 | Bool_t AliTrackResidualsFast::Update() | |
202 | { | |
203 | // Find the alignment parameters | |
204 | // by using the already accumulated | |
205 | // sums | |
206 | TMatrixDSym smatrix(6); | |
207 | TMatrixD sums(1,6); | |
208 | ||
3bbfb291 | 209 | smatrix(0,0) = fSum[0]; |
210 | smatrix(0,1) = smatrix(1,0) = fSum[1]; | |
211 | smatrix(0,2) = smatrix(2,0) = fSum[2]; | |
212 | smatrix(0,3) = smatrix(3,0) = fSum[3]; | |
213 | smatrix(0,4) = smatrix(4,0) = fSum[4]; | |
214 | smatrix(0,5) = smatrix(5,0) = fSum[5]; | |
215 | smatrix(1,1) = fSum[6]; | |
216 | smatrix(1,2) = smatrix(2,1) = fSum[7]; | |
217 | smatrix(1,3) = smatrix(3,1) = fSum[8]; | |
218 | smatrix(1,4) = smatrix(4,1) = fSum[9]; | |
219 | smatrix(1,5) = smatrix(5,1) = fSum[10]; | |
220 | smatrix(2,2) = fSum[11]; | |
221 | smatrix(2,3) = smatrix(3,2) = fSum[12]; | |
222 | smatrix(2,4) = smatrix(4,2) = fSum[13]; | |
223 | smatrix(2,5) = smatrix(5,2) = fSum[14]; | |
224 | smatrix(3,3) = fSum[15]; | |
225 | smatrix(3,4) = smatrix(4,3) = fSum[16]; | |
226 | smatrix(3,5) = smatrix(5,3) = fSum[17]; | |
227 | smatrix(4,4) = fSum[18]; | |
228 | smatrix(4,5) = smatrix(5,4) = fSum[19]; | |
cc345ce3 | 229 | smatrix(5,5) = fSum[20]; |
230 | ||
231 | sums(0,0) = fSum[21]; sums(0,1) = fSum[22]; sums(0,2) = fSum[23]; | |
232 | sums(0,3) = fSum[24]; sums(0,4) = fSum[25]; sums(0,5) = fSum[26]; | |
46ae650f | 233 | |
cc101660 | 234 | |
235 | Int_t fixedparamat[6]={0,0,0,0,0,0}; | |
236 | const Int_t unfixedparam=GetNFreeParam(); | |
237 | Int_t position[6],last=0;//position is of size 6 but only unfiexedparam indeces will be used | |
238 | ||
239 | if(fBFixed[0]==kTRUE){ | |
240 | fixedparamat[0]=1; | |
241 | } | |
242 | else { | |
243 | position[0]=0; | |
244 | last++; | |
245 | } | |
246 | ||
247 | for(Int_t j=1;j<6;j++){ | |
248 | if(fBFixed[j]==kTRUE){ | |
249 | fixedparamat[j]=fixedparamat[j-1]+1; | |
250 | } | |
251 | else { | |
252 | fixedparamat[j]=fixedparamat[j-1]; | |
253 | position[last]=j; | |
254 | last++; | |
255 | } | |
256 | } | |
46ae650f | 257 | |
cc101660 | 258 | TMatrixDSym smatrixRedu(unfixedparam); |
259 | for(Int_t i=0;i<unfixedparam;i++){ | |
260 | for(Int_t j=0;j<unfixedparam;j++){ | |
261 | smatrixRedu(i,j)=smatrix(position[i],position[j]); | |
262 | } | |
263 | } | |
264 | ||
ffc90c05 | 265 | // smatrixRedu.Print(); |
cc101660 | 266 | smatrixRedu.Invert(); |
267 | ||
268 | if (!smatrixRedu.IsValid()) { | |
4c84f475 | 269 | printf("Minimization Failed! \n"); |
270 | return kFALSE; | |
271 | } | |
272 | ||
cc101660 | 273 | TMatrixDSym smatrixUp(6); |
274 | for(Int_t i=0;i<6;i++){ | |
275 | for(Int_t j=0;j<6;j++){ | |
276 | if(fBFixed[i]==kTRUE||fBFixed[j]==kTRUE)smatrixUp(i,j)=0.; | |
277 | else smatrixUp(i,j)=smatrixRedu(i-fixedparamat[i],j-fixedparamat[j]); | |
278 | } | |
279 | } | |
280 | ||
4c84f475 | 281 | Double_t covmatrarray[21]; |
282 | ||
cc101660 | 283 | for(Int_t i=0;i<6;i++){ |
4c84f475 | 284 | for(Int_t j=0;j<=i;j++){ |
cc101660 | 285 | if(fBFixed[i]==kFALSE&&fBFixed[j]==kFALSE){ |
286 | if(TMath::Abs(smatrixUp(i,j)/TMath::Sqrt(TMath::Abs(smatrixUp(i,i)*smatrixUp(j,j))))>1.01)printf("Too large Correlation number!\n"); | |
287 | } | |
288 | covmatrarray[i*(i+1)/2+j]=smatrixUp(i,j); | |
4c84f475 | 289 | } |
cc101660 | 290 | } |
4c84f475 | 291 | |
cc101660 | 292 | TMatrixD res = sums*smatrixUp; |
46ae650f | 293 | fAlignObj->SetPars(res(0,0),res(0,1),res(0,2), |
294 | TMath::RadToDeg()*res(0,3), | |
295 | TMath::RadToDeg()*res(0,4), | |
296 | TMath::RadToDeg()*res(0,5)); | |
4c84f475 | 297 | |
298 | fAlignObj->SetCorrMatrix(covmatrarray); | |
46ae650f | 299 | TMatrixD tmp = res*sums.T(); |
300 | fChi2 = fSumR - tmp(0,0); | |
cc101660 | 301 | fNdf -= unfixedparam; |
4c84f475 | 302 | |
46ae650f | 303 | return kTRUE; |
304 | } | |
4c84f475 | 305 | |
306 | ||
307 |