CommitLineData
55f2eba2 1/**************************************************************************
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
26x(t) = x_{off}(t)+k_N\frac{\Delta{P/T}}{P/T}
27
28where x(t) is the parameter which we observe.
29
30\begin{split}
31x(t)=\frac{\Delta{G}}{G_0} \\
32x(t)=\frac{\Delta{v_d}}{v_{d0}}
33\end{split}
34
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"
26514b4e 52#include "TTreeStream.h"
53#include "TRandom.h"
55f2eba2 54
55
56AliTPCkalmanTime::AliTPCkalmanTime():
57 TNamed(),
58 fState(0),
59 fCovariance(0),
60 fTime(0)
61{
62 //
63 // Default constructor
64 //
65}
66
67AliTPCkalmanTime::AliTPCkalmanTime(Double_t time, Double_t xoff, Double_t k, Double_t sigmaxoff, Double_t sigmak):
68 TNamed(),
69 fState(0),
70 fCovariance(0),
71 fTime(0)
72{
73 //
74 // Default constructor
75 //
76 Init(time,xoff,k,sigmaxoff,sigmak);
77}
78
79
80void AliTPCkalmanTime::Init(Double_t time, Double_t xoff, Double_t k, Double_t sigmaxoff, Double_t sigmak){
81 //
82 // Default constructor
83 //
26514b4e 84 fState = new TMatrixD(2,1);
55f2eba2 85 fCovariance = new TMatrixD(2,2);
86 (*fState)(0,0)= xoff; // offset
87 (*fState)(1,0)= k; // slope of the taylor
88 fTime=time;
89 (*fCovariance)(0,0)=sigmaxoff*sigmaxoff;
90 (*fCovariance)(1,1)=sigmak*sigmak;
91 (*fCovariance)(0,1)=0;
92 (*fCovariance)(1,0)=0;
93}
94
95
26514b4e 96void AliTPCkalmanTime::Propagate(Double_t time, Double_t sigma, TTreeSRedirector *debug){
55f2eba2 97 //
98 // Propagate the Kalman
99 //
100 if (!fCovariance) return;
101 if (!fState) return;
26514b4e 102 Double_t deltaT =time-fTime; //delta time - param2 is the current time
55f2eba2 103 Double_t sigmaT2 =(deltaT*deltaT)*sigma*sigma;
26514b4e 104 if (debug){
105 (*debug)<<"matP"<<
106 "time="<<time<<
107 "fTime="<<fTime<<
108 "sigmaT2="<<sigmaT2<<
109 "cov.="<<fCovariance<<
110 "\n";
111 }
55f2eba2 112 (*fCovariance)(0,0)+=sigmaT2;
113 fTime=time;
114}
115
26514b4e 116void AliTPCkalmanTime::Update(Double_t x, Double_t xerr, Double_t ptratio,TTreeSRedirector *debug){
55f2eba2 117 //
118 //
119 //
120 static TMatrixD matHk(1,2); // vector to mesurement
121 static TMatrixD measR(1,1); // measurement error
122 static TMatrixD matQk(2,2); // prediction noise vector
123 static TMatrixD vecZk(1,1); // measurement
124 //
125 static TMatrixD vecYk(1,1); // Innovation or measurement residual
126 static TMatrixD matHkT(2,1);
127 static TMatrixD matSk(1,1); // Innovation (or residual) covariance
128 static TMatrixD matKk(2,1); // Optimal Kalman gain
129 static TMatrixD mat1(2,2); // update covariance matrix
130 static TMatrixD covXk2(2,2);
131 //
132 //
133 TMatrixD covXk(*fCovariance); // X covariance
134 TMatrixD vecXk(*fState); // X vector
135 //
136 vecZk(0,0) = x; // current mesurement
137 measR(0,0) = xerr*xerr; // measurement error
138 matHk(0,0)=1; matHk(0,1)= ptratio;
139 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
140 //
141 //
142 matHkT=matHk.T(); matHk.T();
143 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
144 matSk.Invert();
145 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
146 vecXk += matKk*vecYk; // updated vector
26514b4e 147 mat1(0,0)=1; mat1(1,1)=1;
148 mat1(1,0)=0; mat1(0,1)=0;
55f2eba2 149 covXk2= (mat1-(matKk*matHk));
55f2eba2 150 //
26514b4e 151 //
152 //
153 covXk = covXk2*covXk;
154 //
155 if (debug){
156 (*debug)<<"matU"<<
157 "time="<<fTime<<
158 "x="<<x<<
159 "xerr="<<xerr<<
160 "pt="<<ptratio<<
161 "matHk.="<<&matHk<<
162 "matHkT.="<<&matHkT<<
163 "matSk.="<<&matSk<<
164 "matKk.="<<&matKk<<
165 "covXk2.="<<&covXk2<<
166 "covXk.="<<&covXk<<
167 "cov.="<<fCovariance<<
168 "vecYk.="<<&vecYk<<
169 "vecXk.="<<&vecXk<<
170 "vec.="<<fState<<
171 "\n";
172 }
55f2eba2 173 (*fCovariance)=covXk;
174 (*fState)=vecXk;
175}
26514b4e 176
177void AliTPCkalmanTime::TestMC(const char * fname){
178 //
179 // Test of the class
180 /*
181 Usage:
182 AliTPCkalmanTime::TestMC("testKalman.root");
183 TFile f("testKalman.root");
184
185 */
186 //
187 TTreeSRedirector *pcstream = new TTreeSRedirector(fname);
188 //
189 const Double_t kp0=0;
190 const Double_t kp1=1;
191 const Int_t klength=10*24*60*60; // 10 days mesurement
192 const Double_t ksigmap0=0.001/(24*60*60.); // 0.005 change in one day
d9382d87 193 const Int_t deltat=5*60; // 5 minutes step
26514b4e 194 const Double_t kmessError=0.0005;
195 AliTPCkalmanTime testKalman;
196
197 for (Int_t iter=0; iter<100;iter++){
198 Double_t sp0 = kp0+gRandom->Gaus(0,0.01); // variable to estimate -offset
199 Double_t sp1 = kp1+gRandom->Gaus(0,0.2); // variable to estimate slope
200 Double_t cp0 = sp0; // variable to estimate
201 Double_t cp1 = sp1;
202
203 //
204 testKalman.Init(0,cp0+gRandom->Gaus(0,0.05),cp1+gRandom->Gaus(0,0.2),0.05,0.2);
205 Double_t dptratio= 0;
206 for (Int_t itime=0; itime<klength; itime+=deltat){
207 dptratio+=gRandom->Gaus(0,0.0005);
208 cp0 +=gRandom->Gaus(0,ksigmap0*deltat);
209 //
210 Double_t vdrift = cp0+dptratio*cp1+gRandom->Gaus(0,kmessError);
211 testKalman.Propagate(itime,ksigmap0,pcstream);
212 Double_t fdrift = (*(testKalman.fState))(0,0) + dptratio*(*(testKalman.fState))(1,0);
213 (*pcstream)<<"drift"<<
214 "iter="<<iter<<
215 "time="<<itime<<
216 "vdrift="<<vdrift<<
217 "fdrift="<<fdrift<<
218 "pt="<<dptratio<<
219 "sp0="<<sp0<<
220 "sp1="<<sp1<<
221 "cp0="<<cp0<<
222 "cp1="<<cp1<<
223 "vecXk.="<<testKalman.fState<<
224 "covXk.="<<testKalman.fCovariance<<
225 "\n";
226 testKalman.Update(vdrift,kmessError,dptratio,pcstream);
227 }
228 }
229 delete pcstream;
230}
231