Adding the AliTPCkalmanTime class (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCkalmanTime.cxx
CommitLineData
55f2eba2 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
20The drift velocity and the gas gain are changing in time.
21The drift velocity and gas gain is a function of many parameters, but not all of
22them are known. We assume that the most important parameters are pressure and temperature
23and the influence of other parameters (gas composition, and electric field) are only
24slowly varying in time and can be expressed by smooth function $x_{off}(t)$:
25\begin{equation}
26x(t) = x_{off}(t)+k_N\frac{\Delta{P/T}}{P/T}
27\end{equation}
28where x(t) is the parameter which we observe.
29\begin{equation}
30\begin{split}
31x(t)=\frac{\Delta{G}}{G_0} \\
32x(t)=\frac{\Delta{v_d}}{v_{d0}}
33\end{split}
34\end{equation}
35
36Kalman 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
42Kalman 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
54AliTPCkalmanTime::AliTPCkalmanTime():
55 TNamed(),
56 fState(0),
57 fCovariance(0),
58 fTime(0)
59{
60 //
61 // Default constructor
62 //
63}
64
65AliTPCkalmanTime::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
78void 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
94void 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
106void 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}