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