STEER/AliTrackResidualsLinear.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
  18 //-----------------------------------------------------------------
  19 //   Implementation of the derived class for track residuals
  20 //   based on linear chi2 minimization
  21 //  The minimization relies on the fact that the alignment parameters
  22 //   (angles and translations) are small.
  23 //   TLinearFitter used for minimization
  24 //   Possibility to fix Paramaters
  25 //   FixParameter()ReleaseParameter();
  26 //   Possibility to define fraction of outliers to be skipped
  27 //
  28 //   marian.ivanov@cern.ch
  29 //
  30 //-----------------------------------------------------------------
  31
  32 #include <TMath.h>
  33 #include <TMinuit.h>
  34 #include <TGeoMatrix.h>
  35
  36 #include "AliLog.h"
  37 #include "AliAlignObj.h"
  38 #include "AliTrackPointArray.h"
  39 #include "AliTrackResidualsLinear.h"
  40 #include "AliAlignObj.h"
  41 #include "TLinearFitter.h"
  42 #include  "TDecompSVD.h"
  43
  44 ClassImp(AliTrackResidualsLinear)
  45
  46 //______________________________________________________________________________
  47 AliTrackResidualsLinear::AliTrackResidualsLinear():
  48   AliTrackResiduals(),
  49   fFitter(0),
  50   fFraction(-1),
  51   fChi2Orig(0)
  52 {
  53   // Default constructor
  54   for (Int_t ipar=0; ipar<6; ipar++){
  55     fBFixed[ipar] = kFALSE;
  56     fFixed[ipar]  = 0;;
  57     fParams[ipar]  = 0;
  58   }
  59   for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;}
  60 }
  61
  62 //______________________________________________________________________________
  63 AliTrackResidualsLinear::AliTrackResidualsLinear(Int_t ntracks):
  64   AliTrackResiduals(ntracks),
  65   fFitter(new TLinearFitter(6,"hyp6")),
  66   fFraction(-1),
  67   fChi2Orig(0)
  68 {
  69   // Constructor
  70   for (Int_t ipar=0; ipar<6; ipar++){
  71     fBFixed[ipar] = kFALSE;
  72     fFixed[ipar]  = 0;
  73     fParams[ipar]  = 0;
  74   }
  75   for (Int_t icov=0; icov<36; icov++){ fCovar[icov]=0;}
  76 }
  77
  78 //______________________________________________________________________________
  79 AliTrackResidualsLinear::AliTrackResidualsLinear(const AliTrackResidualsLinear &res):
  80   AliTrackResiduals(res),
  81   fFitter(new TLinearFitter(*(res.fFitter))),
  82   fFraction(res.fFraction),
  83   fChi2Orig(res.fChi2Orig)
  84 {
  85   // Copy constructor
  86   //..
  87   for (Int_t ipar=0; ipar<6; ipar++){
  88     fBFixed[ipar]  = res.fBFixed[ipar];
  89     fFixed[ipar]   = res.fFixed[ipar];
  90     fParams[ipar]  = res.fParams[ipar];
  91   }
  92   for (Int_t icov=0; icov<36; icov++){ fCovar[icov]= res.fCovar[icov];}
  93   fChi2Orig = res.fChi2Orig;
  94 }
  95
  96 //______________________________________________________________________________
  97 AliTrackResidualsLinear &AliTrackResidualsLinear::operator= (const AliTrackResidualsLinear& res)
  98 {
  99   // Assignment operator
 100  ((AliTrackResiduals *)this)->operator=(res);
 101  return *this;
 102 }
 103 //______________________________________________________________________________
 104 AliTrackResidualsLinear::~AliTrackResidualsLinear()
 105 {
 106   //
 107   //
 108   //
 109   delete fFitter;
 110 }
 111
 112
 113 //______________________________________________________________________________
 114 Bool_t AliTrackResidualsLinear::Minimize()
 115 {
 116   // Implementation of fast linear Chi2 minimizer
 117   // based on TLinear fitter
 118   //
 119   if (!fFitter) fFitter = new TLinearFitter(6,"hyp6");
 120   fFitter->StoreData(kTRUE);
 121   fFitter->ClearPoints();
 122   fChi2Orig = 0;
 123   AliTrackPoint p1,p2;
 124   for (Int_t itrack = 0; itrack < fLast; itrack++) {
 125     if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
 126     for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
 127       fVolArray[itrack]->GetPoint(p1,ipoint);
 128       fTrackArray[itrack]->GetPoint(p2,ipoint);
 129       AddPoints(p1,p2);
 130     }
 131   }
 132   Bool_t isOK = Update();
 133   if (!isOK) return isOK;
 134   //
 135   TGeoHMatrix  matrix;
 136   fAlignObj->GetMatrix(matrix);
 137   return isOK;
 138 }
 139
 140
 141
 142 //______________________________________________________________________________
 143 void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime)
 144 {
 145   //
 146   // add points to linear fitter - option with correlation betwee measurement in different dimensions
 147   // p1      - local point
 148   // pprime  - track extrapolation point
 149   //
 150   Float_t xyz[3],xyzp[3];
 151   Float_t cov[6],covp[6];
 152   p.GetXYZ(xyz,cov); pprime.GetXYZ(xyzp,covp);
 153   //
 154   //
 155   TMatrixD mcov(3,3);   // local point covariance
 156   mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
 157   mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
 158   mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
 159   TMatrixD mcovp(3,3); //  extrapolation point covariance
 160   mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
 161   mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
 162   mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
 163   mcov+=mcovp;
 164   //mcov.Invert();
 165   if (!mcov.IsValid()) return;
 166   TMatrixD mcovBack = mcov;  // for debug purposes
 167   //
 168   // decompose matrix
 169   //
 170   TDecompSVD svd(mcov);              // mcov  = svd.fU * covDiagonal * svd.fV.Invert
 171   if (!svd.Decompose()) return;      // decomposition failed
 172   TMatrixD   matrixV = svd.GetV();   // transformation matrix to diagonalize covariance matrix
 173   Double_t   covDiagonal[3] = {svd.GetSig()[0],svd.GetSig()[1],svd.GetSig()[2]};    // diagonalized covariance matrix
 174   //
 175   // residual vector
 176   TMatrixD  deltaR(3,1);
 177   deltaR(0,0) = (xyzp[0]-xyz[0]);
 178   deltaR(1,0) = (xyzp[1]-xyz[1]);
 179   deltaR(2,0) = (xyzp[2]-xyz[2]);
 180   //
 181   // parametrization matrix
 182   //
 183   TMatrixD        mparam(3,6);
 184   mparam(0,0) = 1;      mparam(1,0) = 0;       mparam(2,0) = 0;            // xshift
 185   mparam(0,1) = 0;      mparam(1,1) = 1;       mparam(2,1) = 0;            // yshift
 186   mparam(0,2) = 0;      mparam(1,2) = 0;       mparam(2,2) = 1;            // zshift
 187   mparam(0,3) = 0;      mparam(1,3) =-xyz[2];  mparam(2,3) = xyz[1];       // x rotation
 188   mparam(0,4) = xyz[2]; mparam(1,4) = 0;       mparam(2,4) =-xyz[0];       // y rotation
 189   mparam(0,5) =-xyz[1]; mparam(1,5) = xyz[0];  mparam(2,5) = 0;            // z rotation
 190   //
 191
 192   TMatrixD  deltaT(matrixV, TMatrixD::kTransposeMult, deltaR);   // tranformed delta
 193   TMatrixD  mparamT(matrixV,TMatrixD::kTransposeMult, mparam);   // tranformed linear transformation
 194   if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>2){
 195     //
 196     // debug part
 197     //
 198     //   covDiag = U^-1 * mcov * V      -- diagonalization of covariance matrix
 199
 200     TMatrixD   matrixU = svd.GetU();                      // transformation matrix to diagonalize covariance matrix
 201     TMatrixD   matrixUI= svd.GetU();
 202     matrixUI.Invert();
 203     //
 204     TMatrixD   test0   = matrixUI*matrixV;                // test matrix - should be unit matrix
 205     TMatrixD   test1   = matrixUI*mcovBack*matrixV;       // test matrix - diagonal - should be diagonal with covDiagonal on diag
 206     TMatrixD   test2   = matrixU.T()*matrixV;             // test ortogonality - shoul be unit
 207     printf("Test matrix 2 - should be unit\n");
 208     test2.Print();
 209     printf("Test matrix 0 - should be unit\n");
 210     test0.Print();
 211     printf("Test matrix 1 - should be diagonal\n");
 212     test1.Print();
 213     printf("Diagonal matrix\n");
 214     svd.GetSig().Print();
 215     printf("Original param matrix\n");
 216     mparam.Print();
 217     printf("Rotated  param matrix\n");
 218     mparamT.Print();
 219     //
 220     //
 221     printf("Trans Matrix:\n");
 222     matrixV.Print();
 223     printf("Delta Orig\n");
 224     deltaR.Print();
 225     printf("Delta Rotated");
 226     deltaT.Print();
 227     //
 228     //
 229   }
 230   //
 231   Double_t sumChi2=0;
 232   for (Int_t idim = 1; idim<3; idim++){
 233     Double_t yf;     // input points to fit in TLinear fitter
 234     Double_t xf[6];  // input points to fit
 235     yf = deltaT(idim,0);
 236     for (Int_t ipar =0; ipar<6; ipar++) xf[ipar] = mparamT(idim,ipar);
 237     if (covDiagonal[idim]>0.){
 238       fFitter->AddPoint(xf,yf, TMath::Sqrt(covDiagonal[idim]));
 239       // accumulate chi2
 240       Double_t chi2       = (yf*yf)/covDiagonal[idim];
 241       fChi2Orig += (yf*yf)/covDiagonal[idim];
 242       if (chi2>100 && AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){
 243         printf("Too big chi2- %f\n",chi2);
 244         printf("Delta Orig\n");
 245         deltaR.Print();
 246         printf("Delta Rotated");
 247         deltaT.Print();
 248         matrixV.Print();
 249         printf("Too big chi2 - End\n");
 250       }
 251       sumChi2+=chi2;
 252     }
 253     else{
 254       printf("Bug\n");
 255     }
 256   }
 257   if (AliLog::GetDebugLevel("","AliTrackResidualsLinear")>1){
 258     TMatrixD matChi0=(mcov.Invert()*deltaR);
 259     TMatrixD matChi2=deltaR.T()*matChi0;
 260     printf("Chi2:\t%f\t%f", matChi2(0,0), sumChi2);
 261   }
 262
 263   fNdf +=2;
 264 }
 265
 266 //______________________________________________________________________________
 267 Bool_t AliTrackResidualsLinear::Update()
 268 {
 269   // Find the alignment parameters
 270   // using TLinear fitter  + fill data containers
 271   //
 272   //
 273   fFitter->Eval();
 274   //
 275   // TLinear fitter put as first parameter offset - fixing parameter shifted by one
 276   //
 277   fFitter->FixParameter(0);
 278   for (Int_t ipar =0; ipar<6; ipar++){
 279     if (fBFixed[ipar])  fFitter->FixParameter(ipar+1,fFixed[ipar]);
 280   }
 281   if (fFraction>0.5) {
 282     fFitter->EvalRobust(fFraction);
 283   }else{
 284     fFitter->Eval();
 285   }
 286   //
 287   fFitter->ReleaseParameter(0);
 288   for (Int_t ipar=0; ipar<6; ipar++) {
 289     if (fBFixed[ipar])  fFitter->ReleaseParameter(ipar+1);
 290   }
 291
 292
 293   //
 294   fChi2 = fFitter->GetChisquare();
 295   fNdf -= 6;
 296   TVectorD vector(7);
 297   fFitter->GetParameters(vector);
 298   fParams[0] = vector[1];
 299   fParams[1] = vector[2];
 300   fParams[2] = vector[3];
 301   fParams[3] = vector[4];
 302   fParams[4] = vector[5];
 303   fParams[5] = vector[6];
 304   TMatrixD covar(7,7);
 305   fFitter->GetCovarianceMatrix(covar);
 306   for (Int_t i0=0; i0 <6; i0++)
 307     for (Int_t j0=0; j0 <6; j0++){
 308       fCovar[i0*6+j0] = covar(i0+1,j0+1);
 309     }
 310   //
 311   fAlignObj->SetPars(fParams[0], fParams[1], fParams[2],
 312                      TMath::RadToDeg()*fParams[3],
 313                      TMath::RadToDeg()*fParams[4],
 314                      TMath::RadToDeg()*fParams[5]);
 315   return kTRUE;
 316 }