STEER/AliTrackFitterStraight.cxx

   1 #include "TMatrixDSym.h"
   2 #include "TMatrixD.h"
   3 #include "AliTrackFitterStraight.h"
   4
   5 ClassImp(AliTrackFitterStraight)
   6
   7 AliTrackFitterStraight::AliTrackFitterStraight():
   8   AliTrackFitter()
   9 {
  10   //
  11   // default constructor
  12   //
  13   fAlpha = 0.;
  14   for (Int_t i=0;i<5;i++) fSumXY[i] = 0;
  15   fSumYY = 0;
  16   for (Int_t i=0;i<5;i++) fSumXZ[i] = 0;
  17   fSumZZ = 0;
  18   fNUsed = 0;
  19   fConv = kFALSE;
  20 }
  21
  22
  23 AliTrackFitterStraight::AliTrackFitterStraight(AliTrackPointArray *array, Bool_t owner):
  24   AliTrackFitter(array,owner)
  25 {
  26   //
  27   // Constructor
  28   //
  29   fAlpha = 0.;
  30   for (Int_t i=0;i<5;i++) fSumXY[i] = 0;
  31   fSumYY = 0;
  32   for (Int_t i=0;i<5;i++) fSumXZ[i] = 0;
  33   fSumZZ = 0;
  34   fNUsed = 0;
  35   fConv = kFALSE;
  36 }
  37
  38 AliTrackFitterStraight::AliTrackFitterStraight(const AliTrackFitterStraight &fitter):
  39   AliTrackFitter(fitter)
  40 {
  41   //
  42   // copy constructor
  43   //
  44   fAlpha = fitter.fAlpha;
  45   for (Int_t i=0;i<5;i++) fSumXY[i]  = fitter.fSumXY[i];
  46   fSumYY = fitter.fSumYY;
  47   for (Int_t i=0;i<5;i++) fSumXZ[i]  = fitter.fSumXZ[i];
  48   fSumZZ = fitter.fSumZZ;
  49   fNUsed = fitter.fNUsed;
  50   fConv = fitter.fConv;
  51 }
  52
  53 //_____________________________________________________________________________
  54 AliTrackFitterStraight &AliTrackFitterStraight::operator =(const AliTrackFitterStraight& fitter)
  55 {
  56   // assignment operator
  57   //
  58   if(this==&fitter) return *this;
  59   ((AliTrackFitter *)this)->operator=(fitter);
  60
  61   fAlpha = fitter.fAlpha;
  62   for (Int_t i=0;i<5;i++) fSumXY[i]  = fitter.fSumXY[i];
  63   fSumYY = fitter.fSumYY;
  64   for (Int_t i=0;i<5;i++) fSumXZ[i]  = fitter.fSumXZ[i];
  65   fSumZZ = fitter.fSumZZ;
  66   fNUsed = fitter.fNUsed;
  67   fConv = fitter.fConv;
  68
  69   return *this;
  70 }
  71
  72 AliTrackFitterStraight::~AliTrackFitterStraight()
  73 {
  74   // destructor
  75   //
  76 }
  77
  78 void AliTrackFitterStraight::Reset()
  79 {
  80   // Reset the track parameters and
  81   // sums
  82   AliTrackFitter::Reset();
  83   fAlpha = 0.;
  84   for (Int_t i=0;i<5;i++) fSumXY[i] = 0;
  85   fSumYY = 0;
  86   for (Int_t i=0;i<5;i++) fSumXZ[i] = 0;
  87   fSumZZ = 0;
  88   fNUsed = 0;
  89   fConv =kFALSE;
  90 }
  91
  92 Bool_t AliTrackFitterStraight::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
  93                                  AliAlignObj::ELayerID layerRangeMin,
  94                                  AliAlignObj::ELayerID layerRangeMax)
  95 {
  96   // Fit the track points. The method takes as an input
  97   // the set of id's (volids) of the volumes in which
  98   // one wants to calculate the residuals.
  99   // The following parameters are used to define the
 100   // range of volumes to be used in the fitting
 101   // As a result two AliTrackPointArray's obects are filled.
 102   // The first one contains the space points with
 103   // volume id's from volids list. The second array of points represents
 104   // the track extrapolations corresponding to the space points
 105   // in the first array. The two arrays can be used to find
 106   // the residuals in the volids and consequently construct a
 107   // chi2 function to be minimized during the alignment
 108   // procedures. For the moment the track extrapolation is taken
 109   // at the space-point reference plane. The reference plane is
 110   // found using the covariance matrix of the point
 111   // (assuming sigma(x)=0 at the reference coordinate system.
 112
 113   Reset();
 114
 115   Int_t npoints = fPoints->GetNPoints();
 116   if (npoints < 2) return kFALSE;
 117
 118   Bool_t isAlphaCalc = kFALSE;
 119   AliTrackPoint p,plocal;
 120
 121   Int_t npVolId = 0;
 122   fNUsed = 0;
 123   Int_t *pindex = new Int_t[npoints];
 124   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 125     {
 126       fPoints->GetPoint(p,ipoint);
 127       UShort_t iVolId = p.GetVolumeID();
 128       if (FindVolId(volIds,iVolId)) {
 129         pindex[npVolId] = ipoint;
 130         npVolId++;
 131       }
 132       if (volIdsFit != 0x0) {
 133         if (!FindVolId(volIdsFit,iVolId)) continue;
 134       }
 135       else {
 136         if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
 137             iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
 138                                                 AliAlignObj::LayerSize(layerRangeMax-
 139                                                                        AliAlignObj::kFirstLayer))) continue;
 140       }
 141       if (!isAlphaCalc) {
 142         fAlpha = p.GetAngle();
 143         isAlphaCalc = kTRUE;
 144       }
 145       plocal = p.Rotate(fAlpha);
 146       AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
 147                TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
 148       fNUsed++;
 149     }
 150
 151   if (fNUsed < 2) {
 152     delete [] pindex;
 153     return kFALSE;
 154   }
 155
 156   Update();
 157
 158   if (!fConv) {
 159     delete [] pindex;
 160     return kFALSE;
 161   }
 162
 163   if (fNUsed < fMinNPoints ) {
 164     delete [] pindex;
 165     return kFALSE;
 166   }
 167
 168   fPVolId = new AliTrackPointArray(npVolId);
 169   fPTrack = new AliTrackPointArray(npVolId);
 170   AliTrackPoint p2;
 171   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
 172     {
 173       Int_t index = pindex[ipoint];
 174       fPoints->GetPoint(p,index);
 175       if (GetPCA(p,p2)) {
 176         Float_t xyz[3],xyz2[3];
 177         p.GetXYZ(xyz); p2.GetXYZ(xyz2);
 178         //      printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);
 179         fPVolId->AddPoint(ipoint,&p);
 180         fPTrack->AddPoint(ipoint,&p2);
 181       }
 182     }
 183
 184   delete [] pindex;
 185
 186   return kTRUE;
 187 }
 188
 189 void AliTrackFitterStraight::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
 190 {
 191   // Straight track fitter
 192   // The method add a point to the sums
 193   // used to extract track parameters
 194
 195   //
 196   // XY part
 197   //
 198   Double_t weight = 1./(sy*sy);
 199   fSumXY[0] +=weight;
 200   fSumXY[1] +=x*weight;      fSumXY[2] +=x*x*weight;
 201   fSumXY[3] +=y*weight;      fSumXY[4] +=x*y*weight;
 202   fSumYY += y*y*weight;
 203   //
 204   // XZ part
 205   //
 206   weight = 1./(sz*sz);
 207   fSumXZ[0] +=weight;
 208   fSumXZ[1] +=x*weight;      fSumXZ[2] +=x*x*weight;
 209   fSumXZ[3] +=z*weight;      fSumXZ[4] +=x*z*weight;
 210   fSumZZ += z*z*weight;
 211 }
 212
 213 void AliTrackFitterStraight::Update(){
 214   //
 215   //  Track fitter update
 216   //
 217   //
 218   for (Int_t i=0;i<6;i++)fParams[i]=0;
 219   fChi2 = 0;
 220   fNdf = 0;
 221   Int_t conv=0;
 222   //
 223   // XY part
 224   //
 225   TMatrixDSym     smatrix(2);
 226   TMatrixD        sums(1,2);
 227   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[2];
 228   smatrix(0,1) = fSumXY[1];
 229   sums(0,0)    = fSumXY[3]; sums(0,1)   =fSumXY[4];
 230   smatrix.Invert();
 231   if (smatrix.IsValid()){
 232     for (Int_t i=0;i<2;i++)
 233       for (Int_t j=0;j<=i;j++){
 234         (*fCov)(i,j)=smatrix(i,j);
 235       }
 236     TMatrixD  res = sums*smatrix;
 237     fParams[0] = res(0,0);
 238     fParams[1] = res(0,1);
 239     TMatrixD  tmp = res*sums.T();
 240     fChi2 += fSumYY - tmp(0,0);
 241     fNdf  += fNUsed - 2;
 242     conv++;
 243   }
 244   //
 245   // XZ part
 246   //
 247   TMatrixDSym     smatrixz(2);
 248   TMatrixD        sumsxz(1,2);
 249   smatrixz(0,0) = fSumXZ[0]; smatrixz(1,1) = fSumXZ[2];
 250   smatrixz(0,1) = fSumXZ[1];
 251   sumsxz(0,0)   = fSumXZ[3]; sumsxz(0,1)   = fSumXZ[4];
 252   smatrixz.Invert();
 253   if (smatrixz.IsValid()){
 254     TMatrixD res = sumsxz*smatrixz;
 255     fParams[2] = res(0,0);
 256     fParams[3] = res(0,1);
 257     fParams[4] = fParams[5] = 0;
 258     for (Int_t i=0;i<2;i++)
 259       for (Int_t j=0;j<=i;j++){
 260         (*fCov)(i+2,j+2)=smatrixz(i,j);
 261       }
 262     TMatrixD  tmp = res*sumsxz.T();
 263     fChi2 += fSumZZ - tmp(0,0);
 264     fNdf  += fNUsed - 2;
 265     conv++;
 266   }
 267
 268   if (conv>1)
 269     fConv =kTRUE;
 270   else
 271     fConv=kFALSE;
 272 }
 273
 274 Double_t AliTrackFitterStraight::GetYat(Double_t x) const {
 275   if (!fConv) return 0.;
 276   return (fParams[0]+x*fParams[1]);
 277 }
 278
 279 Double_t AliTrackFitterStraight::GetDYat(Double_t x) const {
 280   if (!fConv) return 0.;
 281   return fParams[1]+0.*x;
 282 }
 283
 284
 285
 286 Double_t AliTrackFitterStraight::GetZat(Double_t x) const {
 287   if (!fConv) return 0.;
 288   return (fParams[2]+x*fParams[3]);
 289 }
 290
 291 Double_t AliTrackFitterStraight::GetDZat(Double_t x) const {
 292   if (!fConv) return 0.;
 293   return fParams[3]+0.*x;
 294 }
 295
 296 Bool_t AliTrackFitterStraight::GetXYZat(Double_t r, Float_t *xyz) const {
 297   if (!fConv) return kFALSE;
 298   Double_t y = (fParams[0]+r*fParams[1]);
 299   Double_t z = (fParams[2]+r*fParams[3]);
 300
 301   Double_t sin = TMath::Sin(fAlpha);
 302   Double_t cos = TMath::Cos(fAlpha);
 303   xyz[0] = r*cos - y*sin;
 304   xyz[1] = y*cos + r*sin;
 305   xyz[2] = z;
 306
 307   return kTRUE;
 308 }
 309
 310 Bool_t AliTrackFitterStraight::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
 311 {
 312   // Get the closest to a given spacepoint track trajectory point
 313   // Look for details in the description of the Fit() method
 314
 315   if (!fConv) return kFALSE;
 316
 317   // First X and Y coordinates
 318   Double_t sin = TMath::Sin(fAlpha);
 319   Double_t cos = TMath::Cos(fAlpha);
 320   // Track parameters in the global coordinate system
 321   Double_t x0 = -fParams[0]*sin;
 322   Double_t y0 =  fParams[0]*cos;
 323   if ((cos - fParams[1]*sin) == 0) return kFALSE;
 324   Double_t dydx = (fParams[1]*cos + sin)/(cos - fParams[1]*sin);
 325
 326   // Define space-point refence plane
 327   Double_t alphap = p.GetAngle();
 328   Double_t sinp = TMath::Sin(alphap);
 329   Double_t cosp = TMath::Cos(alphap);
 330   Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
 331   //  Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
 332   Double_t x0p= x0*cosp + y0*sinp;
 333   Double_t y0p= y0*cosp - x0*sinp;
 334   if ((cos + dydx*sin) == 0) return kFALSE;
 335   Double_t dydxp = (dydx*cos - sin)/(cos + dydx*sin);
 336   Double_t yprime = y0p + dydxp*(x-x0p);
 337
 338   // Back to the global coordinate system
 339   Double_t xsecond = x*cosp - yprime*sinp;
 340   Double_t ysecond = yprime*cosp + x*sinp;
 341
 342   // Now Z coordinate and track angles
 343   Double_t x2 = xsecond*cos + ysecond*sin;
 344   Double_t zsecond = GetZat(x2);
 345   Double_t dydx2 = fParams[1];
 346   Double_t dzdx = fParams[3];
 347
 348   // Fill the cov matrix of the track extrapolation point
 349   Double_t cov[6] = {0,0,0,0,0,0};
 350   Double_t sigmax = 100*100.;
 351   cov[0] = sigmax;             cov[1] = sigmax*dydx2;      cov[2] = sigmax*dzdx;
 352   cov[3] = sigmax*dydx2*dydx2; cov[4] = sigmax*dydx2*dzdx;
 353   cov[5] = sigmax*dzdx*dzdx;
 354
 355   Float_t  newcov[6];
 356   newcov[0] = cov[0]*cos*cos-
 357             2*cov[1]*sin*cos+
 358               cov[3]*sin*sin;
 359   newcov[1] = cov[1]*(cos*cos-sin*sin)-
 360              (cov[3]-cov[0])*sin*cos;
 361   newcov[2] = cov[2]*cos-
 362               cov[4]*sin;
 363   newcov[3] = cov[0]*sin*sin+
 364             2*cov[1]*sin*cos+
 365               cov[3]*cos*cos;
 366   newcov[4] = cov[4]*cos+
 367               cov[2]*sin;
 368   newcov[5] = cov[5];
 369
 370   p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
 371
 372   return kTRUE;
 373 }