]>
Commit | Line | Data |
---|---|---|
258112da | 1 | |
2 | // See http://en.wikipedia.org/wiki/Kalman_filter | |
3 | // | |
4 | /* | |
5 | .x ~/UliStyle.C | |
6 | gSystem->AddIncludePath("-I$ALICE_ROOT/STAT"); | |
7 | gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); | |
8 | gSystem->AddIncludePath("-I$ALICE_ROOT/TPC"); | |
9 | gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+") | |
10 | gSystem->Load("libSTAT.so"); | |
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 | ||
92 | Double_t cshift = 0; | |
93 | Double_t cPT = 0; | |
94 | Double_t floatS = 0.00005; // | |
95 | Double_t floatTP = 0.0005; // | |
96 | Double_t noiseV = 0.3/250.; // cosmic drift precission | |
97 | Double_t vdrift = 0; | |
98 | ||
99 | AliTPCkalmanTime kalmanTime; | |
100 | ||
101 | ||
102 | void 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 | ||
136 | void 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 | ||
170 | void 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 | } |