8b1e8f29b83b6b4149cb6b037d7e4bcdac831a63
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibAlign.cxx
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>
36 using namespace std;
37
38 ClassImp(AliTPCcalibAlign)
39
40 AliTPCcalibAlign::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
51 AliTPCcalibAlign::~AliTPCcalibAlign() {
52   //
53   // destructor
54   //
55 }
56
57 void 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];
65   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
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
88 void 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];
101   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
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
165 void 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];
178   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
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
242 void 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];
255   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
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
319 void 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
360 Bool_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
389 Bool_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
416 Bool_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 }