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