Update TPCCEda to write output file in parts (to avoid too big files produced in...
[u/mrichter/AliRoot.git] / 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 #include "TTreeStream.h"
53 #include "TRandom.h"
54
55
56 AliTPCkalmanTime::AliTPCkalmanTime():
57   TNamed(),
58   fState(0),
59   fCovariance(0),
60   fTime(0)
61 {
62   //
63   // Default constructor
64   //
65 }
66
67 AliTPCkalmanTime::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
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   fState = new TMatrixD(2,1);
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
96 void AliTPCkalmanTime::Propagate(Double_t time, Double_t sigma, TTreeSRedirector *debug){
97   //
98   // Propagate the Kalman
99   //
100   if (!fCovariance) return;
101   if (!fState) return;
102   Double_t deltaT  =time-fTime;  //delta time - param2 is the current time
103   Double_t sigmaT2 =(deltaT*deltaT)*sigma*sigma; 
104   if (debug){
105     (*debug)<<"matP"<<
106       "time="<<time<<
107       "fTime="<<fTime<<
108       "sigmaT2="<<sigmaT2<<
109       "cov.="<<fCovariance<<
110       "\n";
111   }
112   (*fCovariance)(0,0)+=sigmaT2;
113   fTime=time;  
114 }
115
116 void  AliTPCkalmanTime::Update(Double_t x, Double_t xerr, Double_t ptratio,TTreeSRedirector *debug){
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 
147   mat1(0,0)=1; mat1(1,1)=1; 
148   mat1(1,0)=0; mat1(0,1)=0;
149   covXk2= (mat1-(matKk*matHk));
150   //
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   }
173   (*fCovariance)=covXk;
174   (*fState)=vecXk;
175 }
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   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   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