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 }