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   fSumZZ = 0;
  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   fSumZZ =0;
  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   fSumZZ    = rieman.fSumZZ;
 106   fCapacity = rieman.fN;
 107   fN =rieman.fN;
 108   fCovar = new TMatrixDSym(*(rieman.fCovar));
 109   fConv = rieman.fConv;
 110   fX  = new Double_t[fN];
 111   fY  = new Double_t[fN];
 112   fZ  = new Double_t[fN];
 113   fSy = new Double_t[fN];
 114   fSz = new Double_t[fN];
 115   for (Int_t i=0;i<rieman.fN;i++){
 116     fX[i]   = rieman.fX[i];
 117     fY[i]   = rieman.fY[i];
 118     fZ[i]   = rieman.fZ[i];
 119     fSy[i]  = rieman.fSy[i];
 120     fSz[i]  = rieman.fSz[i];
 121   }
 122 }
 123
 124 AliRieman::~AliRieman()
 125 {
 126   //
 127   // Destructor
 128   //
 129   delete[]fX;
 130   delete[]fY;
 131   delete[]fZ;
 132   delete[]fSy;
 133   delete[]fSz;
 134   delete fCovar;
 135 }
 136
 137 void AliRieman::Reset()
 138 {
 139   //
 140   // Reset all the data members
 141   //
 142   fN=0;
 143   for (Int_t i=0;i<6;i++) fParams[i] = 0;
 144   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
 145   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
 146   fSumZZ =0;
 147   fConv =kFALSE;
 148 }
 149
 150
 151 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
 152 {
 153   //
 154   //  Rieman update
 155   //
 156   //------------------------------------------------------
 157   // XY direction
 158   //
 159   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 160   //
 161   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 162   //
 163   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 164   //
 165   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 166   //
 167   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 168   //
 169   //  final linear equation :   a + u*b +t*c - v =0;
 170   //
 171   // Minimization :
 172   //
 173   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 174   //
 175   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 176   //
 177   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 178   //
 179   if (fN==fCapacity-1) return;  // out of capacity
 180   fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
 181   //
 182   // XY part
 183   //
 184   Double_t  t  =  x*x+y*y;
 185   if (t<2) return;
 186   t            = 1./t;
 187   Double_t  u  =  2.*x*t;
 188   Double_t  v  =  2.*y*t;
 189   Double_t  error = 2.*sy*t;
 190   error *=error;
 191   Double_t weight = 1./error;
 192   fSumXY[0] +=weight;
 193   fSumXY[1] +=u*weight;      fSumXY[2]+=v*weight;  fSumXY[3]+=t*weight;
 194   fSumXY[4] +=u*u*weight;    fSumXY[5]+=t*t*weight;
 195   fSumXY[6] +=u*v*weight;    fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
 196   //
 197   // XZ part
 198   //
 199   weight = 1./(sz*sz);
 200   fSumXZ[0] +=weight;
 201   Double_t r = TMath::Sqrt(x*x+y*y);
 202   fSumXZ[1] +=weight*r;   fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
 203   // Now the auxulary sums
 204   fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
 205   fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
 206   fSumZZ += z*z*weight;
 207   //
 208   fN++;
 209 }
 210
 211
 212
 213
 214
 215 void AliRieman::UpdatePol(){
 216   //
 217   //  Rieman update
 218   //
 219   //
 220   if (fN==0) return;
 221   for (Int_t i=0;i<6;i++)fParams[i]=0;
 222   Int_t conv=0;
 223   //
 224   // XY part
 225   //
 226   TMatrixDSym     smatrix(3);
 227   TMatrixD        sums(1,3);
 228   //
 229   //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
 230   //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
 231   //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
 232
 233   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
 234   smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
 235   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
 236   smatrix.Invert();
 237   if (smatrix.IsValid()){
 238     for (Int_t i=0;i<3;i++)
 239       for (Int_t j=0;j<=i;j++){
 240         (*fCovar)(i,j)=smatrix(i,j);
 241       }
 242     TMatrixD  res = sums*smatrix;
 243     fParams[0] = res(0,0);
 244     fParams[1] = res(0,1);
 245     fParams[2] = res(0,2);
 246     conv++;
 247   }
 248   //
 249   // XZ part
 250   //
 251   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 252   fSumXZ[1] += fSumXZ[5]*Rm1*Rm1;
 253   fSumXZ[2] += fSumXZ[6]*Rm1*Rm1 + fSumXZ[8]*Rm1*Rm1*Rm1*Rm1;
 254   fSumXZ[4] += fSumXZ[7]*Rm1*Rm1;
 255
 256   TMatrixDSym     smatrixz(2);
 257   smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
 258   smatrixz(1,0) = fSumXZ[1];
 259   smatrixz.Invert();
 260   TMatrixD        sumsxz(1,2);
 261   if (smatrixz.IsValid()){
 262     sumsxz(0,0)    = fSumXZ[3];
 263     sumsxz(0,1)    = fSumXZ[4];
 264     TMatrixD res = sumsxz*smatrixz;
 265     fParams[3] = res(0,0);
 266     fParams[4] = res(0,1);
 267     fParams[5] = 0;
 268     for (Int_t i=0;i<2;i++)
 269       for (Int_t j=0;j<=i;j++){
 270         (*fCovar)(i+3,j+3)=smatrixz(i,j);
 271       }
 272     conv++;
 273   }
 274   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 275   //
 276   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 277   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, y = 2*y*t,  D0 = R^2 - x0^2- y0^2
 278   //
 279   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 280   //
 281   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 282   //  final linear equation :   a + u*b +t*c - v =0;
 283   //
 284   //  fParam[0]  = 1/y0
 285   //  fParam[1]  = -x0/y0
 286   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 287   if (conv>1) fConv =kTRUE;
 288   else
 289     fConv=kFALSE;
 290 }
 291
 292 void AliRieman::Update(){
 293   //
 294   //  Rieman update
 295   //
 296   //
 297   if (fN==0) return;
 298   for (Int_t i=0;i<6;i++)fParams[i]=0;
 299   Int_t conv=0;
 300   //
 301   // XY part
 302   //
 303   TMatrixDSym     smatrix(3);
 304   TMatrixD        sums(1,3);
 305   //
 306   //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
 307   //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
 308   //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
 309
 310   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
 311   smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
 312   //
 313   smatrix(1,0) = fSumXY[1];
 314   smatrix(2,0) = fSumXY[3];
 315   smatrix(2,1) = fSumXY[7];
 316
 317   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
 318   //TDecompChol decomp(smatrix);
 319   //decomp.SetTol(1.0e-14);
 320   //smatrix =
 321   //decomp.Invert();
 322   smatrix.Invert();
 323   if (smatrix.IsValid()){
 324     for (Int_t i=0;i<3;i++)
 325       for (Int_t j=0;j<3;j++){
 326         (*fCovar)(i,j)=smatrix(i,j);
 327       }
 328     TMatrixD  res = sums*smatrix;
 329     fParams[0] = res(0,0);
 330     fParams[1] = res(0,1);
 331     fParams[2] = res(0,2);
 332     conv++;
 333   }
 334   if (conv==0){
 335     fConv = kFALSE;   //non converged
 336     return;
 337   }
 338   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
 339     fConv = kFALSE;   //non converged
 340     return;
 341   }
 342   //
 343   // XZ part
 344   //
 345   Double_t x0 = -fParams[1]/fParams[0];
 346   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.);
 347   Double_t sumXZ[9];
 348
 349   for (Int_t i=0;i<9;i++) sumXZ[i]=0;
 350   for (Int_t i=0;i<fN;i++){
 351     Double_t phi  = TMath::ASin((fX[i]-x0)*rm1);
 352     Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 353     Double_t weight = 1/fSz[i];
 354     weight *=weight;
 355     Double_t dphi = (phi-phi0)/rm1;
 356     sumXZ[0] +=weight;
 357     sumXZ[1] +=weight*dphi;
 358     sumXZ[2] +=weight*dphi*dphi;
 359     sumXZ[3] +=weight*fZ[i];
 360     sumXZ[4] +=weight*dphi*fZ[i];
 361
 362   }
 363
 364   TMatrixDSym     smatrixz(2);
 365   TMatrixD        sumsz(1,2);
 366   smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
 367   smatrixz(1,0) = sumXZ[1];
 368   smatrixz.Invert();
 369   if (smatrixz.IsValid()){
 370     sumsz(0,0)    = sumXZ[3];
 371     sumsz(0,1)    = sumXZ[4];
 372     TMatrixD res = sumsz*smatrixz;
 373     fParams[3] = res(0,0);
 374     fParams[4] = res(0,1);
 375     for (Int_t i=0;i<2;i++)
 376       for (Int_t j=0;j<2;j++){
 377         (*fCovar)(i+3,j+3)=smatrixz(i,j);
 378     }
 379     conv++;
 380   }
 381
 382   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 383   //
 384   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 385   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 386   //
 387   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 388   //
 389   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 390   //  final linear equation :   a + u*b +t*c - v =0;
 391   //
 392   //  fParam[0]  = 1/y0
 393   //  fParam[1]  = -x0/y0
 394   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 395   if (conv>1) fConv =kTRUE;
 396   else
 397     fConv=kFALSE;
 398   fChi2Y = CalcChi2Y();
 399   fChi2Z = CalcChi2Z();
 400   fChi2  = fChi2Y +fChi2Z;
 401 }
 402
 403
 404
 405 Double_t AliRieman::GetYat(Double_t x) const {
 406   //
 407   // get y at given x position
 408   //
 409   if (!fConv) return 0.;
 410   Double_t res2 = (x*fParams[0]+fParams[1]);
 411   res2*=res2;
 412   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 413   if (res2<0) return 0;
 414   res2 = TMath::Sqrt(res2);
 415   res2 = (1-res2)/fParams[0];
 416   return res2;
 417 }
 418
 419 Double_t AliRieman::GetDYat(Double_t x) const {
 420   //
 421   // get dy/dx at given x position
 422   //
 423   if (!fConv) return 0.;
 424   Double_t x0 = -fParams[1]/fParams[0];
 425   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
 426   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 427   if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
 428   Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
 429   if (fParams[0]<0) res*=-1.;
 430   return res;
 431 }
 432
 433
 434
 435 Double_t AliRieman::GetZat(Double_t x) const {
 436   //
 437   // get z at given x position
 438   //
 439   if (!fConv) return 0.;
 440   Double_t x0 = -fParams[1]/fParams[0];
 441   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 442   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 443   Double_t phi  = TMath::ASin((x-x0)*rm1);
 444   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 445   Double_t dphi = (phi-phi0);
 446   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
 447   return res;
 448 }
 449
 450 Double_t AliRieman::GetDZat(Double_t x) const {
 451   //
 452   // get dz/dx at given x postion
 453   //
 454   if (!fConv) return 0.;
 455   Double_t x0 = -fParams[1]/fParams[0];
 456   if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 457   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 458   Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
 459   return res;
 460 }
 461
 462
 463 //_____________________________________________________________________________
 464 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
 465 {
 466   //
 467   // Returns position given radius
 468   //
 469   if (!fConv) return kFALSE;
 470   Double_t res2 = (r*fParams[0]+fParams[1]);
 471   res2*=res2;
 472   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 473   if (res2<0) return kFALSE;
 474   res2 = TMath::Sqrt(res2);
 475   res2 = (1-res2)/fParams[0];
 476
 477   if (!fConv) return kFALSE;
 478   Double_t x0 = -fParams[1]/fParams[0];
 479   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 480   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 481   Double_t phi  = TMath::ASin((r-x0)*rm1);
 482   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 483   Double_t dphi = (phi-phi0);
 484   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
 485
 486   Double_t sin = TMath::Sin(alpha);
 487   Double_t cos = TMath::Cos(alpha);
 488   xyz[0] = r*cos - res2*sin;
 489   xyz[1] = res2*cos + r*sin;
 490   xyz[2] = res;
 491
 492   return kTRUE;
 493 }
 494
 495
 496 Double_t AliRieman::GetC() const{
 497   //
 498   // get curvature
 499   //
 500   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 501 }
 502
 503 Double_t AliRieman::CalcChi2Y() const{
 504   //
 505   // calculate chi2 for Y
 506   //
 507   Double_t sumchi2 = 0;
 508   for (Int_t i=0;i<fN;i++){
 509     Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
 510     sumchi2+=chi2*chi2;
 511   }
 512   return sumchi2;
 513 }
 514
 515
 516 Double_t AliRieman::CalcChi2Z() const{
 517   //
 518   // calculate chi2 for Z
 519   //
 520   Double_t sumchi2 = 0;
 521   for (Int_t i=0;i<fN;i++){
 522     Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
 523     sumchi2+=chi2*chi2;
 524   }
 525   return sumchi2;
 526 }
 527
 528 Double_t AliRieman::CalcChi2() const{
 529   //
 530   // sum chi2 in both coord - supposing Y and Z independent
 531   //
 532   return CalcChi2Y()+CalcChi2Z();
 533 }
 534
 535 AliRieman * AliRieman::MakeResiduals() const{
 536   //
 537   // create residual structure - ONLY for debug purposes
 538   //
 539   AliRieman *rieman = new AliRieman(fN);
 540   for (Int_t i=0;i<fN;i++){
 541     rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);
 542   }
 543   return rieman;
 544 }
 545
 546