STEER/AliTrackFitterRieman.cxx

   1 #include "TMatrixDSym.h"
   2 #include "TMatrixD.h"
   3 #include "AliTrackFitterRieman.h"
   4
   5 ClassImp(AliTrackFitterRieman)
   6
   7 AliTrackFitterRieman::AliTrackFitterRieman():
   8   AliTrackFitter()
   9 {
  10   //
  11   // default constructor
  12   //
  13   fAlpha = 0.;
  14   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
  15   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
  16   fConv = kFALSE;
  17 }
  18
  19
  20 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
  21   AliTrackFitter(array,owner)
  22 {
  23   //
  24   // Constructor
  25   //
  26   fAlpha = 0.;
  27   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
  28   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
  29   fConv = kFALSE;
  30 }
  31
  32 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
  33   AliTrackFitter(rieman)
  34 {
  35   //
  36   // copy constructor
  37   //
  38   fAlpha = rieman.fAlpha;
  39   for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
  40   for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
  41   fConv = rieman.fConv;
  42 }
  43
  44 //_____________________________________________________________________________
  45 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
  46 {
  47   // assignment operator
  48   //
  49   if(this==&rieman) return *this;
  50   ((AliTrackFitter *)this)->operator=(rieman);
  51
  52   fAlpha = rieman.fAlpha;
  53   for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
  54   for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
  55   fConv = rieman.fConv;
  56
  57   return *this;
  58 }
  59
  60 AliTrackFitterRieman::~AliTrackFitterRieman()
  61 {
  62   // destructor
  63   //
  64 }
  65
  66 void AliTrackFitterRieman::Reset()
  67 {
  68   // Reset the track parameters and
  69   // rieman sums
  70   AliTrackFitter::Reset();
  71   fAlpha = 0.;
  72   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
  73   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
  74   fConv =kFALSE;
  75 }
  76
  77 Bool_t AliTrackFitterRieman::Fit(UShort_t volId,
  78                                  AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
  79                                  AliAlignObj::ELayerID layerRangeMin,
  80                                  AliAlignObj::ELayerID layerRangeMax)
  81 {
  82   // Fit the track points. The method takes as an input
  83   // the id (volid) of the volume to be skipped from fitting.
  84   // The following two parameters are used to define the
  85   // range of volumes to be used in the fitting
  86   // As a result two AliTrackPointArray's obects are filled.
  87   // The first one contains the space points with
  88   // volume id = volid. The second array of points represents
  89   // the track extrapolation corresponding to the space points
  90   // in the first array. The two arrays can be used to find
  91   // the residuals in the volid and consequently construct a
  92   // chi2 function to be minimized during the alignment
  93   // procedures. For the moment the track extrapolation is taken
  94   // as follows: in XY plane - at the CDA between track circle
  95   // and the space point; in Z - the track extrapolation on the Z
  96   // plane defined by the space point.
  97
  98   pVolId = pTrack = 0x0;
  99   fConv = kFALSE;
 100
 101   Int_t npoints = fPoints->GetNPoints();
 102   if (npoints < 3) return kFALSE;
 103
 104   AliTrackPoint p;
 105   fPoints->GetPoint(p,0);
 106   fAlpha = TMath::ATan2(p.GetY(),p.GetX());
 107   Double_t sin = TMath::Sin(fAlpha);
 108   Double_t cos = TMath::Cos(fAlpha);
 109
 110   Int_t npVolId = 0;
 111   Int_t npused = 0;
 112   Int_t *pindex = new Int_t[npoints];
 113   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 114     {
 115       fPoints->GetPoint(p,ipoint);
 116       UShort_t iVolId = p.GetVolumeID();
 117       if (iVolId == volId) {
 118         pindex[npVolId] = ipoint;
 119         npVolId++;
 120       }
 121       if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
 122           iVolId >= AliAlignObj::LayerToVolUID(layerRangeMax,0)) continue;
 123       Float_t x = p.GetX()*cos + p.GetY()*sin;
 124       Float_t y = p.GetY()*cos - p.GetX()*sin;
 125       AddPoint(x,y,p.GetZ(),1,1);
 126       npused++;
 127     }
 128
 129   if (npused < 3) {
 130     delete [] pindex;
 131     return kFALSE;
 132   }
 133
 134   Update();
 135
 136   if (!fConv) {
 137     delete [] pindex;
 138     return kFALSE;
 139   }
 140
 141   if ((fParams[0] == 0) ||
 142       ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
 143     delete [] pindex;
 144     return kFALSE;
 145   }
 146
 147   pVolId = new AliTrackPointArray(npVolId);
 148   pTrack = new AliTrackPointArray(npVolId);
 149   AliTrackPoint p2;
 150   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
 151     {
 152       Int_t index = pindex[ipoint];
 153       fPoints->GetPoint(p,index);
 154       if (GetPCA(p,p2)) {
 155         pVolId->AddPoint(ipoint,&p);
 156         pTrack->AddPoint(ipoint,&p2);
 157       }
 158     }
 159
 160   delete [] pindex;
 161
 162   return kTRUE;
 163 }
 164
 165 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
 166 {
 167   //
 168   //  Rieman update
 169   //
 170   //------------------------------------------------------
 171   // XY direction
 172   //
 173   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 174   //
 175   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 176   //
 177   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, y = 2*y*t,  D0 = R^2 - x0^2- y0^2
 178   //
 179   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 180   //
 181   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 182   //
 183   //  final linear equation :   a + u*b +t*c - v =0;
 184   //
 185   // Minimization :
 186   //
 187   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 188   //
 189   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 190   //
 191   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 192   //
 193   //
 194   // XY part
 195   //
 196   Double_t  t  =  x*x+y*y;
 197   if (t<2) return;
 198   t            = 1./t;
 199   Double_t  u  =  2.*x*t;
 200   Double_t  v  =  2.*y*t;
 201   Double_t  error = 2.*sy*t;
 202   error *=error;
 203   Double_t weight = 1./error;
 204   fSumXY[0] +=weight;
 205   fSumXY[1] +=u*weight;      fSumXY[2]+=v*weight;  fSumXY[3]+=t*weight;
 206   fSumXY[4] +=u*u*weight;    fSumXY[5]+=t*t*weight;
 207   fSumXY[6] +=u*v*weight;    fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
 208   //
 209   // XZ part
 210   //
 211   weight = 1./sz;
 212   fSumXZ[0] +=weight;
 213   fSumXZ[1] +=weight*x;   fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
 214   fSumXZ[5] +=weight*z;   fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
 215 }
 216
 217 void AliTrackFitterRieman::Update(){
 218   //
 219   //  Rieman update
 220   //
 221   //
 222   for (Int_t i=0;i<6;i++)fParams[i]=0;
 223   Int_t conv=0;
 224   //
 225   // XY part
 226   //
 227   TMatrixDSym     smatrix(3);
 228   TMatrixD        sums(1,3);
 229   //
 230   //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
 231   //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
 232   //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
 233
 234   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
 235   smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
 236   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
 237   smatrix.Invert();
 238   if (smatrix.IsValid()){
 239     for (Int_t i=0;i<3;i++)
 240       for (Int_t j=0;j<=i;j++){
 241         (*fCov)(i,j)=smatrix(i,j);
 242       }
 243     TMatrixD  res = sums*smatrix;
 244     fParams[0] = res(0,0);
 245     fParams[1] = res(0,1);
 246     fParams[2] = res(0,2);
 247     conv++;
 248   }
 249   //
 250   // XZ part
 251   //
 252   TMatrixDSym     smatrixz(3);
 253   smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2];
 254   smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3];
 255   smatrixz(2,2) = fSumXZ[4];
 256   smatrixz.Invert();
 257   if (smatrixz.IsValid()){
 258     sums(0,0)    = fSumXZ[5];
 259     sums(0,1)    = fSumXZ[6];
 260     sums(0,2)    = fSumXZ[7];
 261     TMatrixD res = sums*smatrixz;
 262     fParams[3] = res(0,0);
 263     fParams[4] = res(0,1);
 264     fParams[5] = res(0,2);
 265     for (Int_t i=0;i<3;i++)
 266       for (Int_t j=0;j<=i;j++){
 267         (*fCov)(i+2,j+2)=smatrixz(i,j);
 268     }
 269     conv++;
 270   }
 271
 272   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 273   //
 274   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 275   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, y = 2*y*t,  D0 = R^2 - x0^2- y0^2
 276   //
 277   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 278   //
 279   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 280   //  final linear equation :   a + u*b +t*c - v =0;
 281   //
 282   //  fParam[0]  = 1/y0
 283   //  fParam[1]  = -x0/y0
 284   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 285   if (conv>1) fConv =kTRUE;
 286   else
 287     fConv=kFALSE;
 288 }
 289
 290 Double_t AliTrackFitterRieman::GetYat(Double_t x){
 291   if (!fConv) return 0.;
 292   Double_t res2 = (x*fParams[0]+fParams[1]);
 293   res2*=res2;
 294   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 295   if (res2<0) return 0;
 296   res2 = TMath::Sqrt(res2);
 297   res2 = (1-res2)/fParams[0];
 298   return res2;
 299 }
 300
 301 Double_t AliTrackFitterRieman::GetDYat(Double_t x){
 302   if (!fConv) return 0.;
 303   Double_t x0 = -fParams[1]/fParams[0];
 304   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
 305   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 306   if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0;
 307   Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
 308   if (fParams[0]<0) res*=-1.;
 309   return res;
 310 }
 311
 312
 313
 314 Double_t AliTrackFitterRieman::GetZat(Double_t x){
 315   if (!fConv) return 0.;
 316   Double_t x0 = -fParams[1]/fParams[0];
 317   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 318   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 319   Double_t phi  = TMath::ASin((x-x0)*Rm1);
 320   Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
 321   Double_t dphi = (phi-phi0);
 322   Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
 323   return res;
 324 }
 325
 326 Double_t AliTrackFitterRieman::GetDZat(Double_t x){
 327   if (!fConv) return 0.;
 328   Double_t x0 = -fParams[1]/fParams[0];
 329   if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 330   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 331   Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1));
 332   return res;
 333 }
 334
 335
 336 Double_t AliTrackFitterRieman::GetC(){
 337   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 338 }
 339
 340 Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz){
 341   if (!fConv) return kFALSE;
 342   Double_t res2 = (r*fParams[0]+fParams[1]);
 343   res2*=res2;
 344   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 345   if (res2<0) return kFALSE;
 346   res2 = TMath::Sqrt(res2);
 347   res2 = (1-res2)/fParams[0];
 348
 349   if (!fConv) return kFALSE;
 350   Double_t x0 = -fParams[1]/fParams[0];
 351   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 352   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 353   Double_t phi  = TMath::ASin((r-x0)*Rm1);
 354   Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
 355   Double_t dphi = (phi-phi0);
 356   Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
 357
 358   Double_t sin = TMath::Sin(fAlpha);
 359   Double_t cos = TMath::Cos(fAlpha);
 360   xyz[0] = r*cos - res2*sin;
 361   xyz[1] = res2*cos + r*sin;
 362   xyz[2] = res;
 363
 364   return kTRUE;
 365 }
 366
 367 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
 368 {
 369   // Get the closest to a given spacepoint track trajectory point
 370   // Look for details in the description of the Fit() method
 371
 372   if (!fConv) return kFALSE;
 373
 374   // First X and Y coordinates
 375   Double_t sin = TMath::Sin(fAlpha);
 376   Double_t cos = TMath::Cos(fAlpha);
 377   //  fParam[0]  = 1/y0
 378   //  fParam[1]  = -x0/y0
 379   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 380   if (fParams[0] == 0) return kFALSE;
 381   Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
 382   Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
 383   if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
 384   Double_t R  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
 385                 fParams[0];
 386
 387   Double_t x  = p.GetX();
 388   Double_t y  = p.GetY();
 389   Double_t dR = TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));
 390   Double_t xprime = TMath::Abs(R)*(x-x0)/dR + x0;
 391   Double_t yprime = TMath::Abs(R)*(y-y0)/dR + y0;
 392
 393   // Now Z coordinate
 394   Double_t phi  = TMath::ASin((x-x0)/R);
 395   Double_t phi0 = TMath::ASin((0.-x0)/R);
 396   Double_t dphi = (phi-phi0);
 397   Double_t zprime = fParams[3]+fParams[4]*dphi*R;
 398
 399   p2.SetXYZ(xprime,yprime,zprime);
 400
 401   return kTRUE;
 402 }