New class for ITS coordiante transformations used by AliITSgeom nearly
[u/mrichter/AliRoot.git] / ITS / AliITSAlignmentTrack.cxx
CommitLineData
e8189707 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
27ClassImp(AliITSAlignmentTrack)
28
29//______________________________________________________________________
30AliITSAlignmentTrack::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//______________________________________________________________________
40AliITSAlignmentTrack::~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//______________________________________________________________________
51void 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//______________________________________________________________________
65void 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//______________________________________________________________________
79void 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//______________________________________________________________________
97void 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//______________________________________________________________________
112Double_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//______________________________________________________________________
156Int_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//______________________________________________________________________
192Int_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//______________________________________________________________________
238Int_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//______________________________________________________________________
301void 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//______________________________________________________________________
366void 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}