]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/DriftKalman.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / TPC / CalibMacros / DriftKalman.C
CommitLineData
258112da 1
2// See http://en.wikipedia.org/wiki/Kalman_filter
3//
4/*
5.x ~/UliStyle.C
6gSystem->AddIncludePath("-I$ALICE_ROOT/STAT");
7gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
8gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
9gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
4070f709 10gSystem->Load("libSTAT");
258112da 11
12.L $ALICE_ROOT/TPC/CalibMacros/DriftKalman.C++
13 //
14 //
15 //
16 TFile fgoofie("goofieTree.root");
17 TTree *chainTime = (TTree*)fgoofie.Get("drift");
18 chainTime->SetAlias("mcut","abs(areaF/areaN-1.117)<0.1");
19 chainTime->SetAlias("atpcor","(1000*press/(3.329*(tempF)))-1");
20 chainTime->SetAlias("dvdrift","(((vdrift/2.636)-1))");
21 chainTime->SetAlias("trCorr","0");
22
23 Int_t npoints= chainTime->Draw("time:(dvdrift)-trCorr:(atpcor)","mcut");
24 Double_t *timeV = new Double_t[npoints];
25 Double_t *vdriftV = new Double_t[npoints];
26 Double_t *cPTV = new Double_t[npoints];
27 memcpy(timeV,chainTime->GetV1(),npoints*sizeof(Double_t));
28 memcpy(vdriftV,chainTime->GetV2(),npoints*sizeof(Double_t));
29 memcpy(cPTV,chainTime->GetV3(),npoints*sizeof(Double_t));
30
31 Double_t rms=0.006, corr;
32 covXk(0,0)=0.01; covXk(1,1)=0.1; covXk(2,2)=0.1; covXk(1,0)=0; covXk(0,1)=0;
33 vecXk(0,0)=0; vecXk(1,0)=1; vecXk(2,0)=0; //initial stat
34
35 recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,10/(250*24*60*60.),kFALSE);
36 recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,10/(250*24*60*60.),kTRUE);
37 TFile f("realvdrift.root");
38 drift->Draw("vecXk.fElements[1]:time","isOK","")
39
40
41 //
42 // Real data test usage
43 //
44
45 AliXRDPROOFtoolkit tool;
46 TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
47 chainTime->Lookup();
48 //
49
50 Double_t rms, corr;
51 chainTime->SetAlias("atpcorP","press0/(3.329*(temp0+273.15))-1");
52 chainTime->SetAlias("mcutP","press0>0");
53 chainTime->SetAlias("ncutP","press0<1");
54 chainTime->SetAlias("atpcorG","temp1.fElements[14]/(0.003378*(temp0+273.15))-1");
55 chainTime->SetAlias("mcutG","temp1.fElements[14]>0");
56
57 chainTime->SetAlias("atpcor","(mcutP*atpcorP+(ncutP)*atpcorG)");
58 chainTime->SetAlias("mcut","temp1.fElements[14]>0");
59 chainTime->SetAlias("dvdrift","(fDz/500)");
60
61 fitLinear(chainTime, rms,corr);
62
63
64 Int_t npoints= chainTime->Draw("time:(dvdrift)-trCorr:(atpcor)","mcut&&abs(fDz-500*driftG)<1");
65 Double_t *timeV = new Double_t[npoints];
66 Double_t *vdriftV = new Double_t[npoints];
67 Double_t *cPTV = new Double_t[npoints];
68 memcpy(timeV,chainTime->GetV1(),npoints*sizeof(Double_t));
69 memcpy(vdriftV,chainTime->GetV2(),npoints*sizeof(Double_t));
70 memcpy(cPTV,chainTime->GetV3(),npoints*sizeof(Double_t));
71
72
73 kalmanTime.Init(0,0,-1,0.01,0.01);
74
75 recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,1/(250*24*60*60.),kFALSE);
76 recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,1/(250*24*60*60.),kTRUE);
77
78 TFile f("realvdrift.root");
79 drift->Draw("k.fState.fElements[0]:time","isOK","");
80*/
81
82#include "TMatrixD.h"
83#include "TRandom.h"
84#include "TTreeStream.h"
85#include "TMath.h"
86#include "TString.h"
87#include "TChain.h"
88#include "TStatToolkit.h"
89#include "AliTPCkalmanTime.h"
90
91
92Double_t cshift = 0;
93Double_t cPT = 0;
94Double_t floatS = 0.00005; //
95Double_t floatTP = 0.0005; //
96Double_t noiseV = 0.3/250.; // cosmic drift precission
97Double_t vdrift = 0;
98
99AliTPCkalmanTime kalmanTime;
100
101
102void recVDrift(Int_t npoints, Double_t *timeV, Double_t *vdriftV, Double_t * cPTV, Double_t noise, Double_t sshift){
103 //
104 // chainTime->Draw("time:(dvdrift):(atpcor)","mcut&&abs(((dvdrift))-(atpcor))<0.004");
105 //
106 TTreeSRedirector *pcstream = new TTreeSRedirector("realvdrift.root");
107 //
108 Int_t nused=0;
109 for (Int_t ipoint=0; ipoint<npoints;ipoint++){
110 //
111 if (ipoint==0) kalmanTime.fTime=timeV[ipoint];
112 TMatrixD &vecXk = *(kalmanTime.fState);
113 TMatrixD &covXk = *(kalmanTime.fCovariance);
114 Double_t fdrift = vecXk(0,0)+ cPTV[ipoint]*vecXk(1,0);
115 Double_t ddrift = fdrift-vdriftV[ipoint];
116 Double_t sdrift = TMath::Sqrt(covXk(0,0)+noise*noise);
117 Bool_t isOK = nused<10||TMath::Abs(ddrift)<5*sdrift;
118 kalmanTime.Propagate(timeV[ipoint],sshift);
119 kalmanTime.Update(vdriftV[ipoint],noise,cPTV[ipoint],pcstream);
120 nused++;
121 (*pcstream)<<"drift"<<
122 "ipoint="<<ipoint<<
123 "time="<<timeV[ipoint]<<
124 "isOK="<<isOK<<
125 "noise="<<noise<<
126 "cPT="<<cPTV[ipoint]<<
127 "vdrift="<<vdriftV[ipoint]<<
128 "fdrift="<<fdrift<<
129 "k.="<<&kalmanTime<<
130 "\n";
131 }
132 delete pcstream;
133}
134
135
136void recVDriftSort(Int_t npoints, Double_t *timeV, Double_t *vdriftV, Double_t * cPTV, Double_t noise, Double_t sshift, Bool_t sort){
137 //
138 // sort entries before fit
139 //
140 Int_t *index = new Int_t[npoints];
141 Double_t *tmpsort = new Double_t[npoints];
142 TMath::Sort(npoints, timeV, index,sort);
143 //
144 for (Int_t i=0; i<npoints;i++){
145 tmpsort[i]= timeV[index[i]];
146 }
147 for (Int_t i=0; i<npoints;i++){
148 timeV[i] =tmpsort[i];
149 }
150 //
151 for (Int_t i=0; i<npoints;i++){
152 tmpsort[i]= vdriftV[index[i]];
153 }
154 for (Int_t i=0; i<npoints;i++){
155 vdriftV[i] =tmpsort[i];
156 }
157 //
158 for (Int_t i=0; i<npoints;i++){
159 tmpsort[i]= cPTV[index[i]];
160 }
161 for (Int_t i=0; i<npoints;i++){
162 cPTV[i] =tmpsort[i];
163 }
164 delete [] index;
165 delete [] tmpsort;
166 recVDrift(npoints,timeV,vdriftV,cPTV,noise,sshift);
167}
168
169
170void fitLinear(TChain * chainTime, Double_t& rms, Double_t & corr){
171 //
172 // Fit globaly - make aliases;
173 //
174 TStatToolkit toolkit;
175 Double_t chi2=0;
176 Int_t npoints=0;
177 TVectorD fitParam;
178 TMatrixD covMatrix;
179 TString fitStr="";
180 fitStr+="(atpcor)++";
181 fitStr+="(trigger==1)++";
182 fitStr+="(trigger==2)++";
183 fitStr+="(trigger==4)++";
184 fitStr+="(trigger==8)++";
185
186 TString *strvdpt1 = 0;
187 strvdpt1 = toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut" , chi2,npoints,fitParam,covMatrix);
188 chainTime->SetAlias("driftG",strvdpt1->Data());
189
190 strvdpt1 =toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut&&abs(fDz-500*driftG)<2" , chi2,npoints,fitParam,covMatrix);
191 chainTime->SetAlias("driftG",strvdpt1->Data());
192
193
194 rms = TMath::Sqrt(chi2/npoints);
195 corr = fitParam[1]; // cPT parameter
196 //
197 //
198 fitStr="";
199 fitStr+="(trigger==1)++";
200 fitStr+="(trigger==2)++";
201 fitStr+="(trigger==4)++";
202 fitStr+="(trigger==8)++";
203 TString fitExpr = "(dvdrift)-";
204 fitExpr+=fitParam[1];
205 fitExpr+="* (atpcor)";
206 TString *strCorr = toolkit.FitPlane(chainTime,fitExpr.Data(),fitStr.Data(),"abs(fDz-500*driftG)<1" , chi2,npoints,fitParam,covMatrix);
207 chainTime->SetAlias("trCorr",strCorr->Data());
208
209
210 fitStr="";
211 fitStr+="(atpcor)++";
212 fitStr+="(atpcor)^2++";
213 fitStr+="(trigger==1)++";
214 fitStr+="(trigger==2)++";
215 fitStr+="(trigger==4)++";
216 fitStr+="(trigger==8)++";
217 TString *strvdpt2 = 0;
218 strvdpt2 =toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut&&abs(fDz-500*driftG)<2" , chi2,npoints,fitParam,covMatrix);
219 chainTime->SetAlias("driftG2",strvdpt2->Data());
220
221
222}