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