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