Release version of ITS code
[u/mrichter/AliRoot.git] / ITS / AliITSAlignmentTrack.cxx
1 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2  * See cxx source for full Copyright notice                               */
3 /* $Id$ */
4 /* $Author$ */
5 /* $Date$ */
6 /* $Name$ */
7 /* $Header$ */
8 /*
9    $Log$
10    Revision 1.1.2.2  2000/06/04 16:35:37  nilsen
11    One more try to fix log comments.
12
13    Revision 1.1.2.1  2000/03/02 20:13:52  nilsen
14    A new class useful for ITS detector alignment studdies.
15  */
16 /* $Revision$ */
17
18 // Standard Root Libraries
19 #include <TMath.h>
20
21 // ITS libraries
22 #include "AliITSgeom.h"
23 #include "AliITSAlignmentTrack.h"
24 #include "AliITSstatistics.h"
25 #include "AliITSstatistics2.h"
26
27 ClassImp(AliITSAlignmentTrack)
28
29 //______________________________________________________________________
30 AliITSAlignmentTrack::AliITSAlignmentTrack(){
31
32     ftrack=fnclust=0,fnclustMax=0;
33     ffunc=-1;
34     fclust=0;
35     for(Int_t i=0;i<10;i++) fpar[i]=0.0;
36     fpx=fpy=fpz=fp=fpt=0.0;
37     fChi2=-1.0;
38 }
39 //______________________________________________________________________
40 AliITSAlignmentTrack::~AliITSAlignmentTrack(){
41
42     ftrack=fnclust=0,fnclustMax=0;
43     ffunc=-1;
44     delete[] fclust;
45     fclust=0;
46     for(Int_t i=0;i<10;i++) fpar[i]=0.0;
47     fpx=fpy=fpz=fp=fpt=0.0;
48     fChi2=-1.0;
49 }
50 //______________________________________________________________________
51 void AliITSAlignmentTrack::func0(Double_t *go,Double_t *gi){
52     Double_t x,y,z;
53
54     x = gi[0];
55     y = gi[1];
56     z = gi[2];
57     x = fpar[0]+fpar[1]*z;
58     y = fpar[2]+fpar[3]*z;
59     go[0] = x;
60     go[1] = y;
61     go[2] = z;
62     return;
63 }
64 //______________________________________________________________________
65 void AliITSAlignmentTrack::func1(Double_t *go,Double_t *gi){
66     Double_t x,y,z;
67
68     x = gi[0];
69     y = gi[1];
70     z = gi[2];
71     x = fpar[0]+fpar[1]*y;
72     z = fpar[2]+fpar[3]*y;
73     go[0] = x;
74     go[1] = y;
75     go[2] = z;
76     return;
77 }
78 //______________________________________________________________________
79 void AliITSAlignmentTrack::func2(Double_t *go,Double_t *gi){
80     Double_t x,y,z,r,th;
81
82     x = gi[0];
83     y = gi[1];
84     z = gi[2];
85     th = TMath::ATan2(y-fpar[1],x-fpar[0]);
86     r  = TMath::Hypot(x-fpar[0],y-fpar[1]);
87     if(th<0.0) th += 2.0*TMath::Pi();
88     x = fpar[0]+fpar[2]*TMath::Cos(th);
89     y = fpar[1]+fpar[2]*TMath::Sin(th);
90     z = fpar[3]+fpar[4]*r;
91     go[0] = x;
92     go[1] = y;
93     go[2] = z;
94     return;
95 }
96 //______________________________________________________________________
97 void AliITSAlignmentTrack::func(Double_t *go,Double_t *gi){
98
99     switch (ffunc){
100     case 0:
101         func0(go,gi);
102         return;
103     case 1:
104         func1(go,gi);
105         return;
106     case 2:
107         func2(go,gi);
108         return;
109     } // end switch
110 }
111 //______________________________________________________________________
112 Double_t AliITSAlignmentTrack::ComputeChi2(){
113     Int_t    i,j,k,l;
114     Double_t chi2=0.0,go[3],gi[3];
115
116     switch (ffunc) {
117     case -1:
118         return -1.0;
119         break;
120     case 0:
121         for(i=0;i<fnclust;i++){
122             for(j=0;j<3;j++)gi[j] = fclust[i].fxg[j];
123             func0(go,gi);
124             for(k=0;k<2;k++) for(l=0;l<2;l++) {
125                 //if(k==2 || l==2) continue;
126                 chi2 += (go[k] - gi[k])*fclust[i].fExg[k][l]*(go[l] - gi[l]);
127             } // end for k,l
128         } // end for i
129         break;
130     case 1:
131         for(i=0;i<fnclust;i++){
132             for(j=0;j<3;j++)gi[j] = fclust[i].fxg[j];
133             func1(go,gi);
134             for(k=0;k<3;k++) for(l=0;l<3;l++){
135                 if(k==1 || l==1) continue;
136                 chi2 += (go[k] - gi[k])*fclust[i].fExg[k][l]*(go[l] - gi[l]);
137             } // end for k,l
138         } // end for i
139         break;
140     case 2:
141         for(i=0;i<fnclust;i++){
142             for(j=0;j<3;j++)gi[j] = fclust[i].fxg[j];
143             func0(go,gi);
144             for(k=0;k<3;k++) for(l=0;l<3;l++){
145                 chi2 += (go[k] - gi[k])*fclust[i].fExg[k][l]*(go[l] - gi[l]);
146             } // end for k,l
147         } // end for i
148         break;
149     } // end switch
150
151     fChi2 = (Float_t) chi2;
152
153     return  chi2;
154 }
155 //______________________________________________________________________
156 Int_t AliITSAlignmentTrack::FindCircleCenter(Double_t *xc,Double_t *x1,
157                                             Double_t  *x2,Double_t *x3){
158 ////////////////////////////////////////////////////////////////////////
159 //     This was derived as folows. Given three non-linear points find
160 // the circle that is therefor defined by those three non-linar points.
161 // Assume that the circle is centers at xc,yc and has a radous R. Then
162 // (1) R^2 = (x1-xc)^2 + (y1-yc)^2
163 // (2) R^2 = (x2-xc)^2 + (y2-yc)^2
164 // (3) R^2 = (x3-xc)^2 + (y3-yc)^2.
165 // Now consider the two equations derived from the above
166 // (1) - (2) = x1^2 - x2^2 -2xc(x1-x2) + y1^2 - y2y2 -2yc(y1-y2) = 0
167 // (1) - (3) = x1^2 - x3^2 -2xc(x1-x3) + y1^2 - y3y2 -2yc(y1-y3) = 0
168 // solving these two equations for x0 and y0 gives
169 // xc = +{(y1-y2)*(y1-y3)*(y1-y3)+x1*x1*(y1-y3)+x2*x2*(y3-y1)+x3*x3*(y1-y2)}/2d
170 // yc = -{(x1-x2)*(x1-x3)*(x1-x3)+y1*y1*(x1-x3)+y2*y2*(x3-x1)+y3*y3*(x1-x2)}/2d
171 // with d = (x1-x2)*(y1-y3) - (x1-x3)*(y1-y2)
172 ////////////////////////////////////////////////////////////////////////
173     Double_t d;
174
175     d = (x1[0]-x2[0])*(x1[1]-x3[1]) - (x1[0]-x3[0])*(x1[1]-x2[1]);
176     if(d==0.0) return 0;  // fits to a line!
177
178     xc[0] = (x1[1]-x2[1])*(x1[1]-x3[1])*(x1[1]-x3[1])+
179              x1[0]*x1[0]*(x1[1]-x3[1])+
180              x2[0]*x2[0]*(x3[1]-x1[1])+
181              x3[0]*x3[0]*(x1[1]-x2[1]);
182     xc[1] = (x1[0]-x2[0])*(x1[0]-x3[0])*(x1[0]-x3[0])+
183              x1[1]*x1[1]*(x1[0]-x3[0])+
184              x2[1]*x2[1]*(x3[0]-x1[0])+
185              x3[1]*x3[1]*(x1[0]-x2[0]);
186     xc[0] = +0.5*xc[0]/d;
187     xc[1] = -0.5*xc[1]/d;
188
189     return 1;
190 }
191 //______________________________________________________________________
192 Int_t AliITSAlignmentTrack::FitTrackToLineG(){
193 // Xg = fpar[0] + fpar[1] * Zg
194 // Yg = fpar[2] + fpar[3] * Zg;
195 // Local Variables
196     Int_t   i;
197     Double_t x,y,z,wx,wy;
198     Double_t a,b,c,d;
199     AliITSstatistics2 *sx = new AliITSstatistics2(4);
200     AliITSstatistics2 *sy = new AliITSstatistics2(4);
201     Double_t b0,d0;
202
203     fChi2 = -1.0;
204     ffunc = 0;
205     if(fnclust<3) return -1;
206
207     b  = (fclust[0].fxg[0]-fclust[fnclust].fxg[0])/
208          (fclust[0].fxg[2]-fclust[fnclust].fxg[2]);
209     d  = (fclust[0].fxg[1]-fclust[fnclust].fxg[1])/
210          (fclust[0].fxg[2]-fclust[fnclust].fxg[2]);
211     do{
212         b0 = b;
213         d0 = d;
214         for(i=0;i<fnclust;i++){
215             sx->Reset();
216             sy->Reset();
217             x    = fclust[i].fxg[0];
218             y    = fclust[i].fxg[1];
219             z    = fclust[i].fxg[2];
220             wx   = 1./(1./fclust[i].fExg[0][0] +
221                        (1./fclust[i].fExg[2][2])*b0*b0);// 1.0/rms^2
222             wy   = 1./(1./fclust[i].fExg[1][1] +
223                        (1./fclust[i].fExg[2][2])*d0*d0);// 1.0/rms^2
224             sx->AddValue(x,z,wx);
225             sy->AddValue(y,z,wy);
226        } // end for i
227         fChi2  = sx->FitToLine(a,b);
228         fChi2 += sy->FitToLine(c,d);
229         //} while(fabs(b0-b)<1.E-5 && fabs(d0-d)<1.E-5);
230     } while(TMath::Abs(b0-b)<1.E-5 && TMath::Abs(d0-d)<1.E-5);
231     fpar[0] = a;
232     fpar[1] = b;
233     fpar[2] = c;
234     fpar[3] = d;
235     return 0;
236 }
237 //______________________________________________________________________
238 Int_t AliITSAlignmentTrack::FitTrackToLineL(AliITSgeom *gm){
239 // X = fpar[0] + fpar[1] * y;
240 // Z = fpar[2] + fpar[3] * y;
241 // in the local coordinate system of the detector fclust[0].findex.
242    // Local Variables
243    Int_t   i,j,k;
244    Double_t wx/*,wy*/,wz,Exll[3][3];
245    Double_t a,b,c,d;
246    Double_t xg[3],xl[3],x2g[3],x2l[3];
247    AliITSstatistics2 *Fx  = new AliITSstatistics2(2);
248    AliITSstatistics2 *Fz  = new AliITSstatistics2(2);
249
250    fChi2 = -1.0;
251    ffunc = 1;
252    if(fnclust<3) return -1;
253
254    Int_t Npts = fnclust;
255    for(i=0;i<Npts;i++){
256        for(j=0;j<3;j++)for(k=0;k<3;k++) Exll[j][k] = fclust[i].fExl[j][k];
257        for(j=0;j<3;j++){x2l[j] = fclust[i].fxl[j];}
258        gm->LtoL(fclust[i].findex,fclust[0].findex,x2l,xl);
259        gm->LtoLErrorMatrix(fclust[i].findex,fclust[0].findex,
260                            (Double_t **) fclust[i].fExl,(Double_t **) Exll);
261        wx = Exll[0][0];
262        wz = Exll[2][2];
263        Fx->AddValue(xl[0],xl[1],wx);
264        Fz->AddValue(xl[2],xl[1],wz);
265    } // end for i
266    fChi2  = Fx->FitToLine(a,b);
267    fChi2 += Fz->FitToLine(c,d);
268    fpar[0] = a;
269    fpar[1] = b;
270    fpar[2] = c;
271    fpar[3] = d;
272    // convert to global if posible.
273    xl[0]  = a;
274    xl[1]  = 0.0;
275    xl[2]  = c;
276    x2l[0] = a+b;
277    x2l[1] = 1.0;
278    x2l[2] = c+d;
279    gm->LtoG(fclust[0].findex,xl,xg);
280    gm->LtoG(fclust[0].findex,x2l,x2g);
281    c = xg[2] - x2g[2];
282    if(c!=0.0){
283        b = (xg[0] - x2g[0])/c;
284        d = (xg[1] - x2g[1])/c;
285        a = xg[0] - b*xg[2];
286        c = xg[1] - d*xg[2];
287        fpar[4] = a;
288        fpar[5] = b;
289        fpar[6] = c;
290        fpar[7] = d;
291    }else{
292        fpar[4] = 0.0;
293        fpar[5] = 0.0;
294        fpar[6] = 0.0;
295        fpar[7] = 0.0;
296        return -1;
297    }// end if c!=0.0
298    return 0;
299 }
300 //______________________________________________________________________
301 void AliITSAlignmentTrack::FitToFunction(Int_t n,AliITSgeom *gm){
302     Int_t i,j,k,l;
303     Double_t r,w,xc[3],x1[3],x2[3],x3[3];
304     AliITSstatistics2 *sa = new AliITSstatistics2(4);
305     AliITSstatistics2 *sb = new AliITSstatistics2(4);
306
307     ffunc = -1;
308     if(fnclust<3) return;
309
310     switch (n){
311     case -1:
312         return;
313         break;
314     case 0:
315         ffunc = 0;
316         FitTrackToLineG();
317         return;
318         break;
319     case 1:
320         ffunc = 1;
321         FitTrackToLineL(gm);
322         return;
323         break;
324     case 2:
325         ffunc = 2;
326         sa->Reset();
327         sb->Reset();
328         for(i=0;i<fnclust;i++){
329             for(l=0;l<3;l++) x1[l] = (Double_t) fclust[i].fxg[l];
330             r = TMath::Hypot((Double_t)x1[0],(Double_t)x1[1]);
331             w = (Double_t) (fclust[i].fExg[2][2]);
332             sb->AddValue((Double_t)fclust[i].fxg[2],r,w);
333             if(i<fnclust-2)for(j=i+1;j<fnclust-1;j++)for(k=j+1;k<fnclust;k++){
334                 for(l=0;l<3;l++) x2[l] = (Double_t) fclust[j].fxg[l];
335                 for(l=0;l<3;l++) x3[l] = (Double_t) fclust[k].fxg[l];
336                 w = 0.0;
337                 w += TMath::Hypot((Double_t)fclust[i].fExg[0][0],
338                                   (Double_t)fclust[i].fExg[1][1]);
339                 w += TMath::Hypot(w,
340                            TMath::Hypot((Double_t)fclust[j].fExg[0][0],
341                                         (Double_t)fclust[j].fExg[1][1]));
342                 w += TMath::Hypot(w,
343                            TMath::Hypot((Double_t)fclust[k].fExg[0][0],
344                                         (Double_t)fclust[k].fExg[1][1]));
345                 if(FindCircleCenter(xc,x1,x2,x3)==0){ // Can't find circle
346                     FitToFunction(1,gm); // fit to line.
347                     return;
348                 } // end if
349                 sa->AddValue(xc[1],xc[0],w);
350             } // end for j,k
351         } // end for i
352         fpar[0] = sa->GetMeanX();
353         fpar[1] = sa->GetMeanY();
354         fpar[2] = sb->GetMeanX();
355         fChi2   = sb->FitToLine(fpar[3],fpar[4]);
356         fChi2   = (Float_t) ComputeChi2();
357         return;
358         break;
359     default:
360         return;
361         break;
362     } // end switch
363     return;
364 }
365 //______________________________________________________________________
366 void AliITSAlignmentTrack::Streamer(TBuffer &R__b){
367    // Stream an object of class AliITSAlignmentTrack.
368
369    if (R__b.IsReading()) {
370       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
371       TObject::Streamer(R__b);
372       R__b >> ftrack;
373       R__b >> fnclust;
374       R__b >> fnclustMax;
375 //      R__b >> fclust;
376       R__b >> ffunc;
377       R__b.ReadStaticArray(fpar);
378       R__b >> fChi2;
379       R__b >> fpx;
380       R__b >> fpy;
381       R__b >> fpz;
382       R__b >> fp;
383       R__b >> fpt;
384    } else {
385       R__b.WriteVersion(AliITSAlignmentTrack::IsA());
386       TObject::Streamer(R__b);
387       R__b << ftrack;
388       R__b << fnclust;
389       R__b << fnclustMax;
390 //      R__b << fclust;
391       R__b << ffunc;
392       R__b.WriteArray(fpar, 10);
393       R__b << fChi2;
394       R__b << fpx;
395       R__b << fpy;
396       R__b << fpz;
397       R__b << fp;
398       R__b << fpt;
399    }
400 }