]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibAlign.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
(Marian, Magnus)
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibAlign.cxx
CommitLineData
9318a5b4 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///////////////////////////////////////////////////////////////////////////////
18// //
19// Class to make a internal alignemnt of TPC chambers //
20//
21// Different linear tranformation investigated
22// 12 parameters - arbitrary transformation
23// 9 parameters - scaling fixed to 1
24// 6 parameters -
25////
26//// Only the 12 parameter version works so far...
27////
28////
29////
30
31#include "TLinearFitter.h"
32#include "AliTPCcalibAlign.h"
33#include "AliExternalTrackParam.h"
34
35#include <iostream>
36using namespace std;
37
38ClassImp(AliTPCcalibAlign)
39
40AliTPCcalibAlign::AliTPCcalibAlign()
41 :fFitterArray12(72*72),fFitterArray9(72*72),fFitterArray6(72*72)
42{
43 //
44 // Constructor
45 //
46 for (Int_t i=0;i<72*72;++i) {
47 fPoints[i]=0;
48 }
49}
50
51AliTPCcalibAlign::~AliTPCcalibAlign() {
52 //
53 // destructor
54 //
55}
56
57void AliTPCcalibAlign::Process(const AliExternalTrackParam &tp1,
58 const AliExternalTrackParam &tp2,
59 Int_t s1,Int_t s2) {
60 //
61 // Process function to fill fitters
62 //
63 Double_t t1[5],t2[5];
64 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 65 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3]/*, &dzdx2=t2[4]*/;
9318a5b4 66 x1 =tp1.GetX();
67 y1 =tp1.GetY();
68 z1 =tp1.GetZ();
69 Double_t snp1=tp1.GetSnp();
70 dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
71 Double_t tgl1=tp1.GetTgl();
72 // dz/dx = 1/(cos(theta)*cos(phi))
73 dzdx1=1./TMath::Sqrt((1.+tgl1*tgl1)*(1.-snp1*snp1));
74 // x2 =tp1->GetX();
75 y2 =tp2.GetY();
76 z2 =tp2.GetZ();
77 Double_t snp2=tp2.GetSnp();
78 dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
79 Double_t tgl2=tp2.GetTgl();
80 dzdx1=1./TMath::Sqrt((1.+tgl2*tgl2)*(1.-snp2*snp2));
81
82 Process12(t1,t2,GetOrMakeFitter12(s1,s2));
83 Process9(t1,t2,GetOrMakeFitter12(s1,s2));
84 Process6(t1,t2,GetOrMakeFitter12(s1,s2));
85 ++fPoints[72*s1+s2];
86}
87
88void AliTPCcalibAlign::Process12(Double_t *t1,
89 Double_t *t2,
90 TLinearFitter *fitter) {
91 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
92 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
93 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
94 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
95 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
96 //
97 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
98 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
99 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
100 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 101 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 102
103 // TODO:
104 Double_t sy = 1.;
105 Double_t sz = 1.;
106 Double_t sdydx = 1.;
107 Double_t sdzdx = 1.;
108
109 Double_t p[12];
110 Double_t value;
111
112 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
113 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
114 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
115 for (Int_t i=0; i<12;i++) p[i]=0.;
116 p[3+0] = x1; // a10
117 p[3+1] = y1; // a11
118 p[3+2] = z1; // a12
119 p[9+1] = 1.; // a13
120 p[0+1] = y1*dydx2; // a01
121 p[0+2] = z1*dydx2; // a02
122 p[9+0] = dydx2; // a03
123 value = y2;
124 fitter->AddPoint(p,value,sy);
125
126 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
127 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
128 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
129 for (Int_t i=0; i<12;i++) p[i]=0.;
130 p[6+0] = x1; // a20
131 p[6+1] = y1; // a21
132 p[6+2] = z1; // a22
133 p[9+2] = 1.; // a23
134 p[0+1] = y1*dzdx2; // a01
135 p[0+2] = z1*dzdx2; // a02
136 p[9+0] = dzdx2; // a03
137 value = z2;
138 fitter->AddPoint(p,value,sz);
139
140 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
141 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
142 for (Int_t i=0; i<12;i++) p[i]=0.;
143 p[3+0] = 1.; // a10
144 p[3+1] = dydx1; // a11
145 p[3+2] = dzdx1; // a12
146 p[0+0] = -dydx2; // a00
147 p[0+1] = -dydx1*dydx2; // a01
148 p[0+2] = -dzdx1*dydx2; // a02
149 value = 0.;
150 fitter->AddPoint(p,value,sdydx);
151
152 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
153 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
154 for (Int_t i=0; i<12;i++) p[i]=0.;
155 p[6+0] = 1; // a20
156 p[6+1] = dydx1; // a21
157 p[6+2] = dzdx1; // a22
158 p[0+0] = -dzdx2; // a00
159 p[0+1] = -dydx1*dzdx2; // a01
160 p[0+2] = -dzdx1*dzdx2; // a02
161 value = 0.;
162 fitter->AddPoint(p,value,sdzdx);
163}
164
165void AliTPCcalibAlign::Process9(Double_t *t1,
166 Double_t *t2,
167 TLinearFitter *fitter) {
168 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
169 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
170 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
171 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
172 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
173 //
174 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
175 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
176 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
177 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 178 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 179
180 // TODO:
181 Double_t sy = 1.;
182 Double_t sz = 1.;
183 Double_t sdydx = 1.;
184 Double_t sdzdx = 1.;
185
186 Double_t p[12];
187 Double_t value;
188
189 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
190 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
191 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a03)*dydx2
192 for (Int_t i=0; i<12;i++) p[i]=0.;
193 p[3+0] = x1;
194 p[3+1] = y1;
195 p[3+2] = z1;
196 p[9+1] = 1.;
197 p[0+1] = y1*dydx2;
198 p[0+2] = z1*dydx2;
199 p[9+0] = dydx2;
200 value = y2;
201 fitter->AddPoint(p,value,sy);
202
203 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
204 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
205 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a03)*dzdx2;
206 for (Int_t i=0; i<12;i++) p[i]=0.;
207 p[6+0] = x1;
208 p[6+1] = y1;
209 p[6+2] = z1;
210 p[9+2] = 1.;
211 p[0+1] = y1*dzdx2;
212 p[0+2] = z1*dzdx2;
213 p[9+0] = dzdx2;
214 value = z2;
215 fitter->AddPoint(p,value,sz);
216
217 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
218 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
219 for (Int_t i=0; i<12;i++) p[i]=0.;
220 p[3+0] = 1.;
221 p[3+1] = dydx1;
222 p[3+2] = dzdx1;
223 p[0+0] = -dydx2;
224 p[0+1] = -dydx1*dydx2;
225 p[0+2] = -dzdx1*dydx2;
226 value = 0.;
227 fitter->AddPoint(p,value,sdydx);
228
229 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
230 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
231 for (Int_t i=0; i<12;i++) p[i]=0.;
232 p[6+0] = 1;
233 p[6+1] = dydx1;
234 p[6+2] = dzdx1;
235 p[0+0] = -dzdx2;
236 p[0+1] = -dydx1*dzdx2;
237 p[0+2] = -dzdx1*dzdx2;
238 value = 0.;
239 fitter->AddPoint(p,value,sdzdx);
240}
241
242void AliTPCcalibAlign::Process6(Double_t *t1,
243 Double_t *t2,
244 TLinearFitter *fitter) {
245 // x2 = 1 *x1 +-a01*y1 + 0 +a03
246 // y2 = a01*x1 + 1 *y1 + 0 +a13
247 // z2 = a20*x1 + a21*y1 + 1 *z1 +a23
248 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
249 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
250 //
251 // 1 -a01 0 a03 x -p[0] x p[3]
252 // a10 1 0 a13 ==> p[0] x x p[4]
253 // a20 a21 1 a23 p[1] p[2] x p[5]
254 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
cec17745 255 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
9318a5b4 256
257 // TODO:
258 Double_t sy = 1.;
259 Double_t sz = 1.;
260 Double_t sdydx = 1.;
261 Double_t sdzdx = 1.;
262
263 Double_t p[12];
264 Double_t value;
265
266 // x2 = 1 *x1 +-a01*y1 + 0 +a03
267 // y2 = a01*x1 + 1 *y1 + 0 +a13
268 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
269 for (Int_t i=0; i<12;i++) p[i]=0.;
270 p[3+0] = x1;
271 p[3+1] = y1;
272 p[3+2] = z1;
273 p[9+1] = 1.;
274 p[0+1] = y1*dydx2;
275 p[0+2] = z1*dydx2;
276 p[9+0] = dydx2;
277 value = y2;
278 fitter->AddPoint(p,value,sy);
279
280 // x2 = 1 *x1 +-a01*y1 + 0 + a03
281 // z2 = a20*x1 + a21*y1 + 1 *z1 + a23
282 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 +a03)*dzdx2;
283 for (Int_t i=0; i<12;i++) p[i]=0.;
284 p[6+0] = x1;
285 p[6+1] = y1;
286 p[6+2] = z1;
287 p[9+2] = 1.;
288 p[0+1] = y1*dzdx2;
289 p[0+2] = z1*dzdx2;
290 p[9+0] = dzdx2;
291 value = z2;
292 fitter->AddPoint(p,value,sz);
293
294 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
295 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
296 for (Int_t i=0; i<12;i++) p[i]=0.;
297 p[3+0] = 1.;
298 p[3+1] = dydx1;
299 p[3+2] = dzdx1;
300 p[0+0] = -dydx2;
301 p[0+1] = -dydx1*dydx2;
302 p[0+2] = -dzdx1*dydx2;
303 value = 0.;
304 fitter->AddPoint(p,value,sdydx);
305
306 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
307 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
308 for (Int_t i=0; i<12;i++) p[i]=0.;
309 p[6+0] = 1;
310 p[6+1] = dydx1;
311 p[6+2] = dzdx1;
312 p[0+0] = -dzdx2;
313 p[0+1] = -dydx1*dzdx2;
314 p[0+2] = -dzdx1*dzdx2;
315 value = 0.;
316 fitter->AddPoint(p,value,sdzdx);
317}
318
319void AliTPCcalibAlign::Eval() {
320 TLinearFitter *f;
321 for (Int_t s1=0;s1<72;++s1)
322 for (Int_t s2=0;s2<72;++s2)
323 if ((f=GetFitter12(s1,s2))&&fPoints[72*s1+s2]>12) {
324 // cerr<<s1<<","<<s2<<": "<<fPoints[72*s1+s2]<<endl;
325 if (!f->Eval()) {
326 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
327 }
328 }
329 /*
330
331 fitter->Eval();
332 fitter->Eval();
333 chi212 = align->GetChisquare()/(4.*entries);
334
335 TMatrixD mat(13,13);
336 TVectorD par(13);
337 align->GetParameters(par);
338 align->GetCovarianceMatrix(mat);
339
340 //
341 //
342 for (Int_t i=0; i<12;i++){
343 palign12(i)= par(i+1);
344 for (Int_t j=0; j<12;j++){
345 pcovar12(i,j) = mat(i+1,j+1);
346 pcovar12(i,j) *= chi212;
347 }
348 }
349 //
350 for (Int_t i=0; i<12;i++){
351 psigma12(i) = TMath::Sqrt(pcovar12(i,i));
352 palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i));
353 for (Int_t j=0; j<12;j++){
354 pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j));
355 }
356 }
357 */
358}
359
360Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
361 if (!GetFitter12(s1,s2))
362 return false;
363 else {
364 TVectorD p(12);
365 cerr<<"foo"<<endl;
366 GetFitter12(s1,s2)->GetParameters(p);
367 cerr<<"bar"<<endl;
368 a.ResizeTo(4,4);
369 a[0][0]=p[0];
370 a[0][1]=p[1];
371 a[0][2]=p[2];
372 a[0][3]=p[9];
373 a[1][0]=p[3];
374 a[1][1]=p[4];
375 a[1][2]=p[5];
376 a[1][3]=p[10];
377 a[2][0]=p[6];
378 a[2][1]=p[7];
379 a[2][2]=p[8];
380 a[2][3]=p[11];
381 a[3][0]=0.;
382 a[3][1]=0.;
383 a[3][2]=0.;
384 a[3][3]=1.;
385 return true;
386 }
387}
388
389Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
390 if (!GetFitter9(s1,s2))
391 return false;
392 else {
393 TVectorD p(9);
394 GetFitter9(s1,s2)->GetParameters(p);
395 a.ResizeTo(4,4);
396 a[0][0]=p[0];
397 a[0][1]=p[1];
398 a[0][2]=p[2];
399 a[0][3]=p[9];
400 a[1][0]=p[3];
401 a[1][1]=p[4];
402 a[1][2]=p[5];
403 a[1][3]=p[10];
404 a[2][0]=p[6];
405 a[2][1]=p[7];
406 a[2][2]=p[8];
407 a[2][3]=p[11];
408 a[3][0]=0.;
409 a[3][1]=0.;
410 a[3][2]=0.;
411 a[3][3]=1.;
412 return true;
413 }
414}
415
416Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
417 if (!GetFitter6(s1,s2))
418 return false;
419 else {
420 TVectorD p(6);
421 cerr<<"foo"<<endl;
422 GetFitter6(s1,s2)->GetParameters(p);
423 cerr<<"bar"<<endl;
424 a.ResizeTo(4,4);
425 a[0][0]=p[0];
426 a[0][1]=p[1];
427 a[0][2]=p[2];
428 a[0][3]=p[9];
429 a[1][0]=p[3];
430 a[1][1]=p[4];
431 a[1][2]=p[5];
432 a[1][3]=p[10];
433 a[2][0]=p[6];
434 a[2][1]=p[7];
435 a[2][2]=p[8];
436 a[2][3]=p[11];
437 a[3][0]=0.;
438 a[3][1]=0.;
439 a[3][2]=0.;
440 a[3][3]=1.;
441 return true;
442 }
443}