1 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2 * See cxx source for full Copyright notice */
10 Revision 1.1.2.2 2000/06/04 16:35:37 nilsen
11 One more try to fix log comments.
13 Revision 1.1.2.1 2000/03/02 20:13:52 nilsen
14 A new class useful for ITS detector alignment studdies.
18 // Standard Root Libraries
22 #include "AliITSgeom.h"
23 #include "AliITSAlignmentTrack.h"
24 #include "AliITSstatistics.h"
25 #include "AliITSstatistics2.h"
27 ClassImp(AliITSAlignmentTrack)
29 //______________________________________________________________________
30 AliITSAlignmentTrack::AliITSAlignmentTrack(){
32 ftrack=fnclust=0,fnclustMax=0;
35 for(Int_t i=0;i<10;i++) fpar[i]=0.0;
36 fpx=fpy=fpz=fp=fpt=0.0;
39 //______________________________________________________________________
40 AliITSAlignmentTrack::~AliITSAlignmentTrack(){
42 ftrack=fnclust=0,fnclustMax=0;
46 for(Int_t i=0;i<10;i++) fpar[i]=0.0;
47 fpx=fpy=fpz=fp=fpt=0.0;
50 //______________________________________________________________________
51 void AliITSAlignmentTrack::func0(Double_t *go,Double_t *gi){
57 x = fpar[0]+fpar[1]*z;
58 y = fpar[2]+fpar[3]*z;
64 //______________________________________________________________________
65 void AliITSAlignmentTrack::func1(Double_t *go,Double_t *gi){
71 x = fpar[0]+fpar[1]*y;
72 z = fpar[2]+fpar[3]*y;
78 //______________________________________________________________________
79 void AliITSAlignmentTrack::func2(Double_t *go,Double_t *gi){
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;
96 //______________________________________________________________________
97 void AliITSAlignmentTrack::func(Double_t *go,Double_t *gi){
111 //______________________________________________________________________
112 Double_t AliITSAlignmentTrack::ComputeChi2(){
114 Double_t chi2=0.0,go[3],gi[3];
121 for(i=0;i<fnclust;i++){
122 for(j=0;j<3;j++)gi[j] = fclust[i].fxg[j];
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]);
131 for(i=0;i<fnclust;i++){
132 for(j=0;j<3;j++)gi[j] = fclust[i].fxg[j];
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]);
141 for(i=0;i<fnclust;i++){
142 for(j=0;j<3;j++)gi[j] = fclust[i].fxg[j];
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]);
151 fChi2 = (Float_t) chi2;
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 ////////////////////////////////////////////////////////////////////////
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!
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;
191 //______________________________________________________________________
192 Int_t AliITSAlignmentTrack::FitTrackToLineG(){
193 // Xg = fpar[0] + fpar[1] * Zg
194 // Yg = fpar[2] + fpar[3] * Zg;
197 Double_t x,y,z,wx,wy;
199 AliITSstatistics2 *sx = new AliITSstatistics2(4);
200 AliITSstatistics2 *sy = new AliITSstatistics2(4);
205 if(fnclust<3) return -1;
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]);
214 for(i=0;i<fnclust;i++){
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);
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);
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.
244 Double_t wx/*,wy*/,wz,Exll[3][3];
246 Double_t xg[3],xl[3],x2g[3],x2l[3];
247 AliITSstatistics2 *Fx = new AliITSstatistics2(2);
248 AliITSstatistics2 *Fz = new AliITSstatistics2(2);
252 if(fnclust<3) return -1;
254 Int_t Npts = fnclust;
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);
263 Fx->AddValue(xl[0],xl[1],wx);
264 Fz->AddValue(xl[2],xl[1],wz);
266 fChi2 = Fx->FitToLine(a,b);
267 fChi2 += Fz->FitToLine(c,d);
272 // convert to global if posible.
279 gm->LtoG(fclust[0].findex,xl,xg);
280 gm->LtoG(fclust[0].findex,x2l,x2g);
283 b = (xg[0] - x2g[0])/c;
284 d = (xg[1] - x2g[1])/c;
300 //______________________________________________________________________
301 void AliITSAlignmentTrack::FitToFunction(Int_t n,AliITSgeom *gm){
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);
308 if(fnclust<3) return;
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];
337 w += TMath::Hypot((Double_t)fclust[i].fExg[0][0],
338 (Double_t)fclust[i].fExg[1][1]);
340 TMath::Hypot((Double_t)fclust[j].fExg[0][0],
341 (Double_t)fclust[j].fExg[1][1]));
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.
349 sa->AddValue(xc[1],xc[0],w);
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();
365 //______________________________________________________________________
366 void AliITSAlignmentTrack::Streamer(TBuffer &R__b){
367 // Stream an object of class AliITSAlignmentTrack.
369 if (R__b.IsReading()) {
370 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
371 TObject::Streamer(R__b);
377 R__b.ReadStaticArray(fpar);
385 R__b.WriteVersion(AliITSAlignmentTrack::IsA());
386 TObject::Streamer(R__b);
392 R__b.WriteArray(fpar, 10);