TPC/AliTPCkalmanTime.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   TPC Kalman filter implementation for time dependent variables:
  18
  19
  20 The drift velocity and the gas gain are changing in time.
  21 The drift velocity and gas gain is a function of many parameters, but not all of
  22 them are known. We assume that the most important parameters are pressure and temperature
  23 and the influence of other parameters (gas composition, and electric field) are only
  24 slowly varying in time and can be expressed by smooth function $x_{off}(t)$:
  25 \begin{equation}
  26 x(t) = x_{off}(t)+k_N\frac{\Delta{P/T}}{P/T}
  27 \end{equation}
  28 where x(t) is the parameter which we observe.
  29 \begin{equation}
  30 \begin{split}
  31 x(t)=\frac{\Delta{G}}{G_0}      \\
  32 x(t)=\frac{\Delta{v_d}}{v_{d0}}
  33 \end{split}
  34 \end{equation}
  35
  36 Kalman filter parameters are following:
  37 \begin{itemize}
  38 \item State vector  ($x_{off}(t)$, $k_N$) at given time
  39 \item Covariance matrix
  40 \end{itemize}
  41
  42 Kalman filter implent following functions:
  43 \begin{itemize}
  44 \item Prediction - adding covariance element $\sigma_{xoff}$
  45 \item Update state vector with new measurement vector ($x_t,\frac{\Delta{P/T}}{P/T}$)
  46 \end{itemize}
  47
  48
  49 */
  50
  51 #include "AliTPCkalmanTime.h"
  52
  53
  54 AliTPCkalmanTime::AliTPCkalmanTime():
  55   TNamed(),
  56   fState(0),
  57   fCovariance(0),
  58   fTime(0)
  59 {
  60   //
  61   // Default constructor
  62   //
  63 }
  64
  65 AliTPCkalmanTime::AliTPCkalmanTime(Double_t time, Double_t xoff, Double_t k, Double_t sigmaxoff, Double_t sigmak):
  66   TNamed(),
  67   fState(0),
  68   fCovariance(0),
  69   fTime(0)
  70 {
  71   //
  72   // Default constructor
  73   //
  74   Init(time,xoff,k,sigmaxoff,sigmak);
  75 }
  76
  77
  78 void AliTPCkalmanTime::Init(Double_t time, Double_t xoff, Double_t k, Double_t sigmaxoff, Double_t sigmak){
  79   //
  80   // Default constructor
  81   //
  82   fState = new TMatrixD(3,1);
  83   fCovariance = new TMatrixD(2,2);
  84   (*fState)(0,0)= xoff;  // offset
  85   (*fState)(1,0)= k;     // slope of the taylor
  86   fTime=time;
  87   (*fCovariance)(0,0)=sigmaxoff*sigmaxoff;
  88   (*fCovariance)(1,1)=sigmak*sigmak;
  89   (*fCovariance)(0,1)=0;
  90   (*fCovariance)(1,0)=0;
  91 }
  92
  93
  94 void AliTPCkalmanTime::Propagate(Double_t time, Double_t sigma){
  95   //
  96   // Propagate the Kalman
  97   //
  98   if (!fCovariance) return;
  99   if (!fState) return;
 100   Double_t deltaT  =time-(*fState)(2,0);  //delta time - param2 is the current time
 101   Double_t sigmaT2 =(deltaT*deltaT)*sigma*sigma;
 102   (*fCovariance)(0,0)+=sigmaT2;
 103   fTime=time;
 104 }
 105
 106 void  AliTPCkalmanTime::Update(Double_t x, Double_t xerr, Double_t ptratio){
 107   //
 108   //
 109   //
 110   static TMatrixD matHk(1,2);    // vector to mesurement
 111   static TMatrixD measR(1,1);    // measurement error
 112   static TMatrixD matQk(2,2);    // prediction noise vector
 113   static TMatrixD vecZk(1,1);    // measurement
 114   //
 115   static TMatrixD vecYk(1,1);    // Innovation or measurement residual
 116   static TMatrixD matHkT(2,1);
 117   static TMatrixD matSk(1,1);    // Innovation (or residual) covariance
 118   static TMatrixD matKk(2,1);    // Optimal Kalman gain
 119   static TMatrixD mat1(2,2);     // update covariance matrix
 120   static TMatrixD covXk2(2,2);
 121   //
 122   //
 123   TMatrixD covXk(*fCovariance);    // X covariance
 124   TMatrixD vecXk(*fState);         // X vector
 125   //
 126   vecZk(0,0) = x;                // current mesurement
 127   measR(0,0) = xerr*xerr;        // measurement error
 128   matHk(0,0)=1;  matHk(0,1)= ptratio;
 129   vecYk = vecZk-matHk*vecXk;     // Innovation or measurement residual
 130   //
 131   //
 132   matHkT=matHk.T(); matHk.T();
 133   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 134   matSk.Invert();
 135   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 136   vecXk += matKk*vecYk;                    //  updated vector
 137   mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1;
 138   covXk2= (mat1-(matKk*matHk));
 139   covXk =  covXk2*covXk;
 140   //
 141   (*fCovariance)=covXk;
 142   (*fState)=vecXk;
 143 }