STEER/AliRieman.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 /* $Id$ */
  17 // Class for global helix fit of a track
  18 // Author: M.Ivanov
  19 // The method uses the following idea:
  20 //------------------------------------------------------
  21 // XY direction
  22 //
  23 //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
  24 //
  25 //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
  26 //
  27 //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
  28 //
  29 //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
  30 //
  31 //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
  32 //
  33 //  final linear equation :   a + u*b +t*c - v =0;
  34 //
  35 // Minimization :
  36 //
  37 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
  38 //
  39 // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
  40 //
  41 // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
  42
  43
  44 #include "TMatrixDSym.h"
  45 //#include "TDecompChol.h"
  46 #include "TMatrixD.h"
  47 #include "AliRieman.h"
  48
  49 ClassImp(AliRieman)
  50
  51
  52
  53 AliRieman::AliRieman(){
  54   //
  55   // default constructor
  56   //
  57   for (Int_t i=0;i<6;i++) fParams[i] = 0;
  58   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
  59   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
  60
  61   fCapacity = 0;
  62   fN =0;
  63   fX  = 0;
  64   fY  = 0;
  65   fZ  = 0;
  66   fSy = 0;
  67   fSz = 0;
  68   fChi2  = 0;
  69   fChi2Y = 0;
  70   fChi2Z = 0;
  71   fCovar = 0;
  72   fConv = kFALSE;
  73 }
  74
  75
  76 AliRieman::AliRieman(Int_t capacity){
  77   //
  78   // default constructor
  79   //
  80   for (Int_t i=0;i<6;i++) fParams[i] = 0;
  81   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
  82   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
  83
  84   fCapacity = capacity;
  85   fN =0;
  86   fX  = new Double_t[fCapacity];
  87   fY  = new Double_t[fCapacity];
  88   fZ  = new Double_t[fCapacity];
  89   fSy = new Double_t[fCapacity];
  90   fSz = new Double_t[fCapacity];
  91   fCovar = new TMatrixDSym(6);
  92   fChi2  = 0;
  93   fChi2Y = 0;
  94   fChi2Z = 0;
  95   fConv = kFALSE;
  96 }
  97
  98 AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){
  99   //
 100   // copy constructor
 101   //
 102   for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
 103   for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
 104   for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
 105   fCapacity = rieman.fN;
 106   fN =rieman.fN;
 107   fCovar = new TMatrixDSym(*(rieman.fCovar));
 108   fConv = rieman.fConv;
 109   fX  = new Double_t[fN];
 110   fY  = new Double_t[fN];
 111   fZ  = new Double_t[fN];
 112   fSy = new Double_t[fN];
 113   fSz = new Double_t[fN];
 114   for (Int_t i=0;i<rieman.fN;i++){
 115     fX[i]   = rieman.fX[i];
 116     fY[i]   = rieman.fY[i];
 117     fZ[i]   = rieman.fZ[i];
 118     fSy[i]  = rieman.fSy[i];
 119     fSz[i]  = rieman.fSz[i];
 120   }
 121 }
 122
 123 AliRieman::~AliRieman()
 124 {
 125   //
 126   // Destructor
 127   //
 128   delete[]fX;
 129   delete[]fY;
 130   delete[]fZ;
 131   delete[]fSy;
 132   delete[]fSz;
 133   delete fCovar;
 134 }
 135
 136 void AliRieman::Reset()
 137 {
 138   //
 139   // Reset all the data members
 140   //
 141   fN=0;
 142   for (Int_t i=0;i<6;i++) fParams[i] = 0;
 143   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
 144   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
 145   fConv =kFALSE;
 146 }
 147
 148
 149 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
 150 {
 151   //
 152   //  Rieman update
 153   //
 154   //------------------------------------------------------
 155   // XY direction
 156   //
 157   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 158   //
 159   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 160   //
 161   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 162   //
 163   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 164   //
 165   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 166   //
 167   //  final linear equation :   a + u*b +t*c - v =0;
 168   //
 169   // Minimization :
 170   //
 171   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 172   //
 173   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 174   //
 175   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 176   //
 177   if (fN==fCapacity-1) return;  // out of capacity
 178   fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
 179   //
 180   // XY part
 181   //
 182   Double_t  t  =  x*x+y*y;
 183   if (t<2) return;
 184   t            = 1./t;
 185   Double_t  u  =  2.*x*t;
 186   Double_t  v  =  2.*y*t;
 187   Double_t  error = 2.*sy*t;
 188   error *=error;
 189   Double_t weight = 1./error;
 190   fSumXY[0] +=weight;
 191   fSumXY[1] +=u*weight;      fSumXY[2]+=v*weight;  fSumXY[3]+=t*weight;
 192   fSumXY[4] +=u*u*weight;    fSumXY[5]+=t*t*weight;
 193   fSumXY[6] +=u*v*weight;    fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
 194   //
 195   // XZ part
 196   //
 197   weight = 1./sz;
 198   fSumXZ[0] +=weight;
 199   fSumXZ[1] +=weight*x;   fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
 200   fSumXZ[5] +=weight*z;   fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
 201   //
 202   fN++;
 203 }
 204
 205
 206
 207
 208
 209 void AliRieman::UpdatePol(){
 210   //
 211   //  Rieman update
 212   //
 213   //
 214   if (fN==0) return;
 215   for (Int_t i=0;i<6;i++)fParams[i]=0;
 216   Int_t conv=0;
 217   //
 218   // XY part
 219   //
 220   TMatrixDSym     smatrix(3);
 221   TMatrixD        sums(1,3);
 222   //
 223   //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
 224   //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
 225   //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
 226
 227   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
 228   smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
 229   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
 230   smatrix.Invert();
 231   if (smatrix.IsValid()){
 232     for (Int_t i=0;i<3;i++)
 233       for (Int_t j=0;j<=i;j++){
 234         (*fCovar)(i,j)=smatrix(i,j);
 235       }
 236     TMatrixD  res = sums*smatrix;
 237     fParams[0] = res(0,0);
 238     fParams[1] = res(0,1);
 239     fParams[2] = res(0,2);
 240     conv++;
 241   }
 242   //
 243   // XZ part
 244   //
 245   TMatrixDSym     smatrixz(3);
 246   smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2];
 247   smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3];
 248   smatrixz(2,2) = fSumXZ[4];
 249   smatrixz.Invert();
 250   if (smatrixz.IsValid()){
 251     sums(0,0)    = fSumXZ[5];
 252     sums(0,1)    = fSumXZ[6];
 253     sums(0,2)    = fSumXZ[7];
 254     TMatrixD res = sums*smatrixz;
 255     fParams[3] = res(0,0);
 256     fParams[4] = res(0,1);
 257     fParams[5] = res(0,2);
 258     for (Int_t i=0;i<3;i++)
 259       for (Int_t j=0;j<=i;j++){
 260         (*fCovar)(i+2,j+2)=smatrixz(i,j);
 261     }
 262     conv++;
 263   }
 264
 265   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 266   //
 267   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 268   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, y = 2*y*t,  D0 = R^2 - x0^2- y0^2
 269   //
 270   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 271   //
 272   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 273   //  final linear equation :   a + u*b +t*c - v =0;
 274   //
 275   //  fParam[0]  = 1/y0
 276   //  fParam[1]  = -x0/y0
 277   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 278   if (conv>1) fConv =kTRUE;
 279   else
 280     fConv=kFALSE;
 281 }
 282
 283 void AliRieman::Update(){
 284   //
 285   //  Rieman update
 286   //
 287   //
 288   if (fN==0) return;
 289   for (Int_t i=0;i<6;i++)fParams[i]=0;
 290   Int_t conv=0;
 291   //
 292   // XY part
 293   //
 294   TMatrixDSym     smatrix(3);
 295   TMatrixD        sums(1,3);
 296   //
 297   //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
 298   //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
 299   //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
 300
 301   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
 302   smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
 303   //
 304   smatrix(1,0) = fSumXY[1];
 305   smatrix(2,0) = fSumXY[3];
 306   smatrix(2,1) = fSumXY[7];
 307
 308   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
 309   //TDecompChol decomp(smatrix);
 310   //decomp.SetTol(1.0e-14);
 311   //smatrix =
 312   //decomp.Invert();
 313   smatrix.Invert();
 314   if (smatrix.IsValid()){
 315     for (Int_t i=0;i<3;i++)
 316       for (Int_t j=0;j<3;j++){
 317         (*fCovar)(i,j)=smatrix(i,j);
 318       }
 319     TMatrixD  res = sums*smatrix;
 320     fParams[0] = res(0,0);
 321     fParams[1] = res(0,1);
 322     fParams[2] = res(0,2);
 323     conv++;
 324   }
 325   if (conv==0){
 326     fConv = kFALSE;   //non converged
 327     return;
 328   }
 329   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
 330     fConv = kFALSE;   //non converged
 331     return;
 332   }
 333   //
 334   // XZ part
 335   //
 336   Double_t x0 = -fParams[1]/fParams[0];
 337   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.);
 338   Double_t sumXZ[9];
 339
 340   for (Int_t i=0;i<9;i++) sumXZ[i]=0;
 341   for (Int_t i=0;i<fN;i++){
 342     Double_t phi  = TMath::ASin((fX[i]-x0)*rm1);
 343     Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 344     Double_t weight = 1/fSz[i];
 345     weight *=weight;
 346     Double_t dphi = (phi-phi0)/rm1;
 347     sumXZ[0] +=weight;
 348     sumXZ[1] +=weight*dphi;
 349     sumXZ[2] +=weight*dphi*dphi;
 350     sumXZ[3] +=weight*fZ[i];
 351     sumXZ[4] +=weight*dphi*fZ[i];
 352
 353   }
 354
 355   TMatrixDSym     smatrixz(2);
 356   TMatrixD        sumsz(1,2);
 357   smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
 358   smatrixz(1,0) = sumXZ[1];
 359   smatrixz.Invert();
 360   if (smatrixz.IsValid()){
 361     sumsz(0,0)    = sumXZ[3];
 362     sumsz(0,1)    = sumXZ[4];
 363     TMatrixD res = sumsz*smatrixz;
 364     fParams[3] = res(0,0);
 365     fParams[4] = res(0,1);
 366     for (Int_t i=0;i<2;i++)
 367       for (Int_t j=0;j<2;j++){
 368         (*fCovar)(i+3,j+3)=smatrixz(i,j);
 369     }
 370     conv++;
 371   }
 372
 373   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 374   //
 375   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 376   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 377   //
 378   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 379   //
 380   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 381   //  final linear equation :   a + u*b +t*c - v =0;
 382   //
 383   //  fParam[0]  = 1/y0
 384   //  fParam[1]  = -x0/y0
 385   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 386   if (conv>1) fConv =kTRUE;
 387   else
 388     fConv=kFALSE;
 389   fChi2Y = CalcChi2Y();
 390   fChi2Z = CalcChi2Z();
 391   fChi2  = fChi2Y +fChi2Z;
 392 }
 393
 394
 395
 396 Double_t AliRieman::GetYat(Double_t x) const {
 397   //
 398   // get y at given x position
 399   //
 400   if (!fConv) return 0.;
 401   Double_t res2 = (x*fParams[0]+fParams[1]);
 402   res2*=res2;
 403   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 404   if (res2<0) return 0;
 405   res2 = TMath::Sqrt(res2);
 406   res2 = (1-res2)/fParams[0];
 407   return res2;
 408 }
 409
 410 Double_t AliRieman::GetDYat(Double_t x) const {
 411   //
 412   // get dy/dx at given x position
 413   //
 414   if (!fConv) return 0.;
 415   Double_t x0 = -fParams[1]/fParams[0];
 416   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
 417   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 418   if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
 419   Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
 420   if (fParams[0]<0) res*=-1.;
 421   return res;
 422 }
 423
 424
 425
 426 Double_t AliRieman::GetZat(Double_t x) const {
 427   //
 428   // get z at given x position
 429   //
 430   if (!fConv) return 0.;
 431   Double_t x0 = -fParams[1]/fParams[0];
 432   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 433   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 434   Double_t phi  = TMath::ASin((x-x0)*rm1);
 435   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 436   Double_t dphi = (phi-phi0);
 437   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
 438   return res;
 439 }
 440
 441 Double_t AliRieman::GetDZat(Double_t x) const {
 442   //
 443   // get dz/dx at given x postion
 444   //
 445   if (!fConv) return 0.;
 446   Double_t x0 = -fParams[1]/fParams[0];
 447   if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 448   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 449   Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
 450   return res;
 451 }
 452
 453
 454 Double_t AliRieman::GetC() const{
 455   //
 456   // get curvature
 457   //
 458   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 459 }
 460
 461 Double_t AliRieman::CalcChi2Y() const{
 462   //
 463   // calculate chi2 for Y
 464   //
 465   Double_t sumchi2 = 0;
 466   for (Int_t i=0;i<fN;i++){
 467     Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
 468     sumchi2+=chi2*chi2;
 469   }
 470   return sumchi2;
 471 }
 472
 473
 474 Double_t AliRieman::CalcChi2Z() const{
 475   //
 476   // calculate chi2 for Z
 477   //
 478   Double_t sumchi2 = 0;
 479   for (Int_t i=0;i<fN;i++){
 480     Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
 481     sumchi2+=chi2*chi2;
 482   }
 483   return sumchi2;
 484 }
 485
 486 Double_t AliRieman::CalcChi2() const{
 487   //
 488   // sum chi2 in both coord - supposing Y and Z independent
 489   //
 490   return CalcChi2Y()+CalcChi2Z();
 491 }
 492
 493 AliRieman * AliRieman::MakeResiduals() const{
 494   //
 495   // create residual structure - ONLY for debug purposes
 496   //
 497   AliRieman *rieman = new AliRieman(fN);
 498   for (Int_t i=0;i<fN;i++){
 499     rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);
 500   }
 501   return rieman;
 502 }
 503
 504