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