]>
Commit | Line | Data |
---|---|---|
0bfb86a0 | 1 | /* |
2 | marian.ivanov@cern.ch | |
3 | ||
4 | AliTPCkalmanFit: | |
5 | ||
6 | Kalman filter(s) for fitting of the tracks together with calibration/transformation | |
7 | parameters. | |
8 | ||
9 | Correction/Transformation are currently described by set (TObjArray) of primitive | |
10 | correction/transformatio - AliTPCTransformation. Currently we assume that transformation | |
11 | comute (in first order). AliTPCTransformation describe general non linear transformation. | |
12 | ||
13 | Current calibration parameters and covariance stored (fCalibParam, fCalibCovar). | |
14 | ||
15 | Currenly only linear track model implemented. | |
16 | Fits to be implemented: | |
17 | 0. Plane fitting (Laser CE) | |
18 | 1. Primary vertex fitting. | |
19 | 2. Propagation in magnetic field and fit of planes | |
20 | ||
21 | ||
22 | ||
23 | How to use it - see AliTPCkalmanFit::Test function | |
24 | ||
25 | Simple test (see AliTPCkalmanFit::Test) | |
26 | AliTPCTransformation::BuildBasicFormulas(); | |
27 | AliTPCkalmanFit *kalmanFit0 = AliTPCkalmanFit::Test(2000); | |
28 | TFile f("kalmanfitTest.root"); | |
29 | ||
30 | Transformation visualization: | |
31 | Transforamtion can be visualized using TFn (TF1,TF2 ...) and using tree->Draw() | |
32 | e.g: | |
33 | kalmanFit0->SetInstance(kalmanFit0); // | |
34 | kalmanFit0->InitTransformation(); // | |
2497be21 | 35 | // |
7dfa0700 | 36 | TF2 fxRZdz("fxRZdz","sign(y)*1000*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,0,x,0,y)-AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,0,x,0,y-1))",85,245,-250,250); |
0bfb86a0 | 37 | fxRZdz->Draw(""); |
38 | ||
2497be21 | 39 | TF2 fxRZ("fxRZ","sign(y)*10*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y))",85,245,-250,250); |
40 | fxRZ->Draw(""); | |
41 | ||
0bfb86a0 | 42 | |
43 | ||
44 | */ | |
45 | ||
46 | #include "TRandom.h" | |
47 | #include "TMath.h" | |
48 | #include "TBits.h" | |
49 | #include "TFormula.h" | |
50 | #include "TF1.h" | |
51 | #include "TLinearFitter.h" | |
52 | #include "TFile.h" | |
53 | #include "THnSparse.h" | |
54 | #include "TAxis.h" | |
55 | ||
56 | ||
57 | #include "TTreeStream.h" | |
58 | #include "AliTrackPointArray.h" | |
59 | #include "AliLog.h" | |
60 | #include "AliTPCTransformation.h" | |
61 | #include "AliTPCkalmanFit.h" | |
62 | ||
63 | ClassImp(AliTPCkalmanFit) | |
64 | ||
65 | AliTPCkalmanFit* AliTPCkalmanFit::fgInstance = 0; | |
66 | ||
67 | AliTPCkalmanFit::AliTPCkalmanFit(): | |
68 | TNamed(), | |
69 | fCalibration(0), | |
70 | fCalibParam(0), | |
71 | fCalibCovar(0), | |
72 | fLinearParam(0), | |
73 | fLinearCovar(0), | |
74 | fLastTimeStamp(-1), | |
75 | fCurrentAlpha(0), //! current rotation frame | |
76 | fCA(0), //! cosine of current angle | |
77 | fSA(0) //! sinus of current angle | |
78 | { | |
79 | // | |
80 | // Default constructor | |
81 | // | |
c64c82db | 82 | for (Int_t ihis=0; ihis<12; ihis++){ |
0bfb86a0 | 83 | fLinearTrackDelta[ihis]=0; |
84 | fLinearTrackPull[ihis]=0; | |
85 | } | |
86 | } | |
87 | ||
88 | void AliTPCkalmanFit::InitTransformation(){ | |
89 | // | |
90 | // Initialize pointers to the transforamtion functions | |
91 | // | |
92 | // | |
93 | Int_t ncalibs = fCalibration->GetEntries(); | |
94 | for (Int_t icalib=0;icalib<ncalibs; icalib++){ | |
95 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
96 | transform->Init(); | |
97 | } | |
98 | } | |
99 | ||
0e9efcbe | 100 | void AliTPCkalmanFit::Add(const AliTPCkalmanFit * kalman){ |
101 | // | |
102 | // | |
103 | // | |
104 | Update(kalman); | |
c64c82db | 105 | for (Int_t i=0;i<12;i++){ |
0e9efcbe | 106 | if (fLinearTrackDelta[i] && kalman->fLinearTrackDelta[i]){ |
107 | fLinearTrackDelta[i]->Add(kalman->fLinearTrackDelta[i]); | |
108 | } | |
109 | if (fLinearTrackPull[i] && kalman->fLinearTrackPull[i]){ | |
110 | fLinearTrackPull[i]->Add(kalman->fLinearTrackPull[i]); | |
111 | } | |
112 | } | |
113 | ||
114 | } | |
115 | ||
116 | ||
0bfb86a0 | 117 | void AliTPCkalmanFit::Init(){ |
118 | // | |
119 | // Initialize parameter vector and covariance matrix | |
120 | // To be called after initialization of all of the transformations | |
121 | // | |
122 | // | |
123 | Int_t ncalibs = fCalibration->GetEntries(); | |
124 | fCalibParam = new TMatrixD(ncalibs,1); | |
125 | fCalibCovar = new TMatrixD(ncalibs,ncalibs); | |
126 | for (Int_t icalib=0;icalib<ncalibs; icalib++){ | |
127 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
98a4cc77 | 128 | (*fCalibParam)(icalib,0) = transform->GetParam(); |
0bfb86a0 | 129 | for (Int_t jcalib=0;jcalib<ncalibs; jcalib++){ |
130 | if (icalib!=jcalib) (*fCalibCovar)(icalib,jcalib)= 0; | |
98a4cc77 | 131 | if (icalib==jcalib) (*fCalibCovar)(icalib,jcalib) = transform->GetSigma()*transform->GetSigma(); |
0bfb86a0 | 132 | } |
133 | } | |
134 | // | |
135 | // Build QA histograms | |
136 | // | |
137 | Double_t mpi = TMath::Pi(); | |
138 | // | |
139 | // delta alpha y0 z0 ky kz | |
140 | Int_t binsQA[6] = {300, 36*4, 30, 25, 20, 20}; | |
141 | Double_t xminQA[6] = {-0.5, -mpi, -120, -250, -1, -1.}; | |
142 | Double_t xmaxQA[6] = { 0.5, mpi, 120, 250, 1, 1.}; | |
143 | TString axisName[6]={"#Delta", | |
144 | "#alpha", | |
145 | "y_{0}", | |
146 | "z_{0}", | |
147 | "tan(#phi)", | |
148 | "tan(theta)"}; | |
149 | // | |
c64c82db | 150 | TString deltaName[12]={"#Delta_{y}(cm)", |
151 | "100*#Delta_{#phi}(cm)", | |
152 | "100^{2}dy0^{2}/dx0^{2}(cm)", | |
153 | "100^{2}dy1^{2}/dx1^{2}(cm)", | |
154 | "#Delta_{z}(cm)", | |
155 | "100*#Delta_{#theta}(cm)", | |
156 | "100^{2}*dz0^{2}/dx0^{2}(cm)", | |
157 | "100^{2}*dz1^{2}/dx1^{2}(cm)", | |
158 | "RMSy_{0} (cm)", | |
159 | "RMSy_{1} (cm)", | |
160 | "RMSz_{0} (cm)", | |
161 | "RMSz_{1} (cm)"}; | |
162 | ||
163 | TString pullName[12]={"#Delta_{y}(unit)", | |
164 | "100*#Delta_{#phi}(unit)", | |
165 | "100^{2}dy0^{2}/dx0^{2}(unit)", | |
166 | "100^{2}dy1^{2}/dx1^{2}(unit)", | |
167 | "#Delta_{z}(unit)", | |
168 | "100*#Delta_{#theta}(unit)", | |
169 | "100^{2}*dz0^{2}/dx0^{2}(unit)", | |
170 | "100^{2}*dz1^{2}/dx1^{2}(unit)" | |
171 | "RMSy_{0} (cm)", | |
172 | "RMSy_{1} (cm)", | |
173 | "RMSz_{0} (cm)", | |
174 | "RMSz_{1} (cm)"}; | |
175 | // | |
176 | // | |
177 | // | |
178 | for (Int_t ihis=0; ihis<12; ihis++){ | |
0bfb86a0 | 179 | fLinearTrackDelta[ihis]=0; |
180 | fLinearTrackPull[ihis]=0; | |
181 | xminQA[0]=-0.5; xmaxQA[0] = 0.5; | |
182 | fLinearTrackDelta[ihis] = new THnSparseS(deltaName[ihis],deltaName[ihis], 6, binsQA,xminQA, xmaxQA); | |
183 | xminQA[0]=-10; xmaxQA[0] = 10; | |
184 | fLinearTrackPull[ihis] = new THnSparseS(pullName[ihis],pullName[ihis], 6, binsQA,xminQA, xmaxQA); | |
185 | for (Int_t iaxis=1; iaxis<6; iaxis++){ | |
186 | fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]); | |
187 | fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]); | |
188 | fLinearTrackPull[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]); | |
189 | fLinearTrackPull[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]); | |
190 | } | |
191 | fLinearTrackDelta[ihis]->GetAxis(0)->SetName(deltaName[ihis]); | |
192 | fLinearTrackDelta[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]); | |
193 | fLinearTrackPull[ihis]->GetAxis(0)->SetName(deltaName[ihis]); | |
194 | fLinearTrackPull[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]); | |
195 | } | |
196 | ||
197 | ||
198 | } | |
199 | ||
0e9efcbe | 200 | void AliTPCkalmanFit::SetStatus(const char * mask, Bool_t setOn, Bool_t isOr){ |
201 | // | |
202 | // 0. To activate all transforamtion call SetStatus(0,kTRUE) | |
203 | // 1. To disable everything SetStatus(0,kFALSE) | |
204 | // 2. To activate/desactivate SetStatus("xxx",kTRUE/kFALSE,kFALSE) | |
205 | Int_t ncalibs = fCalibration->GetEntries(); | |
206 | if (mask==0) { | |
207 | for (Int_t i=0; i<ncalibs;i++){ | |
208 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i); | |
209 | transform->SetActive(setOn); | |
210 | } | |
211 | } | |
212 | else{ | |
213 | for (Int_t i=0; i<ncalibs;i++){ | |
214 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i); | |
215 | TString strName(transform->GetName()); | |
216 | if (strName.Contains(mask)){ | |
217 | if (!isOr) transform->SetActive( transform->IsActive() && setOn); | |
218 | if (isOr) transform->SetActive( transform->IsActive() || setOn); | |
219 | } | |
220 | } | |
221 | } | |
222 | } | |
223 | ||
224 | ||
225 | void AliTPCkalmanFit::Update(const AliTPCkalmanFit * kalman){ | |
226 | // | |
227 | // Update Kalman filter | |
228 | // | |
229 | Int_t ncalibs = fCalibration->GetEntries(); | |
230 | TMatrixD vecXk=*fCalibParam; // X vector | |
231 | TMatrixD covXk=*fCalibCovar; // X covariance | |
232 | TMatrixD &vecZk = *(kalman->fCalibParam); | |
233 | TMatrixD &measR = *(kalman->fCalibCovar); | |
234 | ||
235 | TMatrixD matHk(ncalibs,ncalibs); // vector to mesurement | |
236 | TMatrixD vecYk(ncalibs,1); // Innovation or measurement residual | |
237 | TMatrixD matHkT(ncalibs,ncalibs); // helper matrix Hk transpose | |
238 | TMatrixD matSk(ncalibs,ncalibs); // Innovation (or residual) covariance | |
239 | TMatrixD matKk(ncalibs,ncalibs); // Optimal Kalman gain | |
240 | TMatrixD covXk2(ncalibs,ncalibs); // helper matrix | |
241 | TMatrixD covXk3(ncalibs,ncalibs); // helper matrix | |
242 | // | |
243 | for (Int_t i=0;i<ncalibs;i++){ | |
244 | for (Int_t j=0;j<ncalibs;j++) matHk(i,j)=0; | |
245 | matHk(i,i)=1; | |
246 | } | |
247 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual | |
248 | matHkT=matHk.T(); matHk.T(); | |
249 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
250 | matSk.Invert(); | |
251 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
252 | vecXk += matKk*vecYk; // updated vector | |
253 | covXk2= (matHk-(matKk*matHk)); | |
254 | covXk3 = covXk2*covXk; | |
255 | covXk = covXk3; | |
256 | (*fCalibParam) = vecXk; | |
257 | (*fCalibCovar) = covXk; | |
258 | } | |
259 | ||
260 | ||
261 | ||
262 | ||
0bfb86a0 | 263 | |
2497be21 | 264 | void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug, Float_t scalingRMSY, Float_t scalingRMSZ){ |
0bfb86a0 | 265 | // |
266 | // | |
267 | ||
268 | if (!fCalibParam) { | |
269 | AliError("Kalman Fit not initialized"); | |
270 | return; | |
271 | } | |
272 | // | |
273 | // 1. Make initial track hypothesis | |
274 | // | |
275 | TLinearFitter lfitY(2,"pol1"); | |
276 | TLinearFitter lfitZ(2,"pol1"); | |
277 | TVectorD vecZ(2); | |
278 | TVectorD vecY(2); | |
279 | // | |
280 | lfitY.StoreData(kTRUE); | |
281 | lfitZ.StoreData(kTRUE); | |
282 | lfitY.ClearPoints(); | |
283 | lfitZ.ClearPoints(); | |
284 | Int_t npoints = points.GetNPoints(); | |
285 | if (npoints<2) return; | |
286 | const Double_t kFac=npoints*npoints*100; | |
287 | const Double_t kOff0=0.01*0.01; | |
288 | const Double_t kOff1=kOff0/(250.*250.); | |
289 | // | |
290 | // 0. Make seed | |
291 | // | |
292 | // AliTrackPointArray pointsc(points); | |
293 | //ApplyCalibration(&pointsc, -1.); | |
294 | // | |
295 | // 1.a Choosing reference rotation alpha | |
296 | // | |
297 | fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]); | |
298 | fCA = TMath::Cos(fCurrentAlpha); | |
299 | fSA = TMath::Sin(fCurrentAlpha); | |
300 | // | |
301 | // 1.b Fit the track in the rotated frame - MakeSeed | |
302 | // | |
c64c82db | 303 | for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){ |
0bfb86a0 | 304 | Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint]; |
305 | Double_t ry = -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint]; | |
306 | Double_t rz = points.GetZ()[ipoint]; | |
307 | lfitY.AddPoint(&rx,ry,1); | |
308 | lfitZ.AddPoint(&rx,rz,1); | |
309 | } | |
310 | lfitY.Eval(); | |
311 | lfitZ.Eval(); | |
312 | lfitY.GetParameters(vecY); | |
313 | lfitZ.GetParameters(vecZ); | |
314 | // | |
315 | // 2. Set initial parameters and the covariance matrix | |
316 | // | |
317 | Int_t ncalibs = fCalibration->GetEntries(); | |
318 | if (!fLinearParam) { | |
319 | fLinearParam = new TMatrixD(ncalibs+4,1); | |
320 | fLinearCovar = new TMatrixD(ncalibs+4,ncalibs+4); | |
321 | } | |
322 | // | |
323 | for (Int_t i=0; i<ncalibs+4;i++){ | |
324 | (*fLinearParam)(i,0)=0; | |
325 | if (i<ncalibs) (*fLinearParam)(i,0) = (*fCalibParam)(i,0); | |
326 | for (Int_t j=0; j<ncalibs+4;j++){ | |
327 | (*fLinearCovar)(i,j)=0; | |
328 | if (i<ncalibs&&j<ncalibs) (*fLinearCovar)(i,j) = (*fCalibCovar)(i,j); | |
329 | } | |
330 | } | |
331 | Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints(); | |
332 | Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints(); | |
333 | (*fLinearParam)(ncalibs+0,0) = lfitY.GetParameter(0); | |
334 | (*fLinearCovar)(ncalibs+0,ncalibs+0)= lfitY.GetCovarianceMatrixElement(0,0)*chi2Y*kFac+kOff0; | |
335 | (*fLinearParam)(ncalibs+1,0) = lfitY.GetParameter(1); | |
336 | (*fLinearCovar)(ncalibs+1,ncalibs+1)= lfitY.GetCovarianceMatrixElement(1,1)*chi2Y*kFac+kOff1; | |
337 | // | |
338 | (*fLinearParam)(ncalibs+2,0) = lfitZ.GetParameter(0); | |
339 | (*fLinearCovar)(ncalibs+2,ncalibs+2)= lfitZ.GetCovarianceMatrixElement(0,0)*chi2Z*kFac+kOff0; | |
340 | (*fLinearParam)(ncalibs+3,0) = lfitZ.GetParameter(1); | |
341 | (*fLinearCovar)(ncalibs+3,ncalibs+3)= lfitZ.GetCovarianceMatrixElement(1,1)*chi2Z*kFac+kOff1; | |
342 | // | |
343 | // Fit thetrack together with correction | |
344 | // | |
345 | AliTrackPoint point; | |
c64c82db | 346 | for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){ |
0bfb86a0 | 347 | // |
348 | if (!points.GetPoint(point,ipoint)) continue; | |
349 | Double_t erry2 = chi2Y; | |
350 | Double_t errz2 = chi2Z; | |
351 | //set error - it is hack | |
352 | Float_t *cov = (Float_t*) point.GetCov(); | |
2497be21 | 353 | cov[1] = (erry2+kOff0)*scalingRMSY; |
354 | cov[2] = (errz2+kOff0)*scalingRMSZ; | |
0bfb86a0 | 355 | UpdateLinear(point,debug); |
356 | if (!points.GetPoint(point,npoints-1-ipoint)) continue; | |
357 | //set error - it is hack | |
358 | cov = (Float_t*) point.GetCov(); | |
2497be21 | 359 | cov[1] = (erry2+kOff0)*scalingRMSY; |
360 | cov[2] = (errz2+kOff0)*scalingRMSZ; | |
0bfb86a0 | 361 | UpdateLinear(point,debug); |
362 | } | |
363 | // | |
364 | // save current param and covariance | |
365 | for (Int_t i=0; i<ncalibs;i++){ | |
366 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i); | |
98a4cc77 | 367 | transform->SetParam( (*fLinearParam)(i,0)); |
0bfb86a0 | 368 | (*fCalibParam)(i,0) = (*fLinearParam)(i,0); |
369 | for (Int_t j=0; j<ncalibs;j++){ | |
370 | (*fCalibCovar)(i,j) = (*fLinearCovar)(i,j); | |
371 | } | |
372 | } | |
373 | if (debug){ // dump debug info | |
374 | (*debug)<<"fitLinear"<< | |
375 | "vecY.="<<&vecY<< | |
376 | "vecZ.="<<&vecZ<< | |
377 | "chi2Y="<<chi2Y<< | |
378 | "chi2Z="<<chi2Z<< | |
379 | "lP.="<<fLinearParam<< | |
380 | "cP.="<<fCalibParam<< | |
381 | "lC.="<<fLinearCovar<< | |
382 | "cC.="<<fCalibCovar<< | |
383 | "\n"; | |
384 | } | |
385 | } | |
386 | ||
387 | void AliTPCkalmanFit::DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){ | |
388 | // | |
389 | // Dump the track parameters before and after current calibration | |
390 | // | |
391 | // Track is divided to two parts - | |
392 | // X mean is defined as middle point | |
393 | // | |
394 | ||
395 | if (!fCalibParam) { | |
396 | AliError("Kalman Fit not initialized"); | |
397 | return; | |
398 | } | |
399 | // | |
400 | // | |
401 | // | |
402 | TLinearFitter *fitters[16]; | |
403 | TVectorD *params[16]; | |
404 | TVectorD *errs[16]; | |
405 | TVectorD chi2N(16); | |
406 | Int_t npoints = points.GetNPoints(); | |
407 | AliTrackPointArray pointsTrans(points); | |
408 | ApplyCalibration(&pointsTrans,-1.); | |
409 | for (Int_t ifit=0; ifit<8;ifit++){ | |
410 | fitters[ifit] = new TLinearFitter(2,"pol1"); | |
411 | params[ifit] = new TVectorD(2); | |
412 | fitters[ifit+8]= new TLinearFitter(3,"pol2"); | |
413 | params[ifit+8] = new TVectorD(3); | |
414 | errs[ifit] = new TVectorD(2); | |
415 | errs[ifit+8] = new TVectorD(3); | |
416 | } | |
417 | // | |
418 | // calculate mean x point and corrdinate frame | |
419 | // | |
420 | fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]); | |
421 | fCA = TMath::Cos(fCurrentAlpha); | |
422 | fSA = TMath::Sin(fCurrentAlpha); | |
423 | Double_t xmean=0, sum=0; | |
424 | for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){ | |
425 | Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint]; | |
426 | xmean+=rx; | |
427 | sum++; | |
428 | } | |
429 | xmean/=sum; | |
430 | // | |
431 | for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){ | |
432 | Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint]; | |
433 | Double_t ry = -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint]; | |
434 | Double_t rz = points.GetZ()[ipoint]; | |
435 | // | |
436 | Double_t rxT = fCA*pointsTrans.GetX()[ipoint]+fSA*pointsTrans.GetY()[ipoint]; | |
437 | Double_t ryT = -fSA*pointsTrans.GetX()[ipoint]+fCA*pointsTrans.GetY()[ipoint]; | |
438 | Double_t rzT = pointsTrans.GetZ()[ipoint]; | |
439 | if (rx>xmean){ | |
440 | fitters[0]->AddPoint(&rxT, ryT,1); | |
441 | fitters[2]->AddPoint(&rxT, rzT,1); | |
442 | fitters[4]->AddPoint(&rx, ry,1); | |
443 | fitters[6]->AddPoint(&rx, rz,1); | |
444 | fitters[8]->AddPoint(&rxT, ryT,1); | |
445 | fitters[10]->AddPoint(&rxT, rzT,1); | |
446 | fitters[12]->AddPoint(&rx, ry,1); | |
447 | fitters[14]->AddPoint(&rx, rz,1); | |
448 | }else{ | |
449 | fitters[1]->AddPoint(&rxT, ryT,1); | |
450 | fitters[3]->AddPoint(&rxT, rzT,1); | |
451 | fitters[5]->AddPoint(&rx, ry,1); | |
452 | fitters[7]->AddPoint(&rx, rz,1); | |
453 | fitters[9]->AddPoint(&rxT, ryT,1); | |
454 | fitters[11]->AddPoint(&rxT, rzT,1); | |
455 | fitters[13]->AddPoint(&rx, ry,1); | |
456 | fitters[15]->AddPoint(&rx, rz,1); | |
457 | } | |
458 | } | |
459 | for (Int_t ifit=0;ifit<16;ifit++){ | |
460 | fitters[ifit]->Eval(); | |
461 | fitters[ifit]->GetParameters(*(params[ifit])); | |
462 | fitters[ifit]->GetErrors(*(errs[ifit])); | |
463 | chi2N[ifit] = TMath::Sqrt(fitters[ifit]->GetChisquare()/(fitters[ifit]->GetNpoints()-1)); | |
464 | (*(errs[ifit]))[0]*=chi2N[ifit]; | |
465 | (*(errs[ifit]))[1]*=chi2N[ifit]; | |
466 | if (ifit>=8) (*(errs[ifit]))[2]*=chi2N[ifit]; // second derivative | |
467 | } | |
468 | if (debug){ | |
469 | (*debug)<<"dumpLinear"<< | |
470 | "alpha="<<fCurrentAlpha<< | |
471 | "xmean="<<xmean<< | |
472 | "y0T.="<<params[0]<< | |
473 | "y1T.="<<params[1]<< | |
474 | "z0T.="<<params[2]<< | |
475 | "z1T.="<<params[3]<< | |
476 | "y0O.="<<params[4]<< | |
477 | "y1O.="<<params[5]<< | |
478 | "z0O.="<<params[6]<< | |
479 | "z1O.="<<params[7]<< | |
480 | "y0T2.="<<params[8]<< | |
481 | "y1T2.="<<params[9]<< | |
482 | "z0T2.="<<params[10]<< | |
483 | "z1T2.="<<params[11]<< | |
484 | "y0O2.="<<params[12]<< | |
485 | "y1O2.="<<params[13]<< | |
486 | "z0O2.="<<params[14]<< | |
487 | "z1O2.="<<params[15]<< | |
488 | "y0TErr.="<<errs[0]<< | |
489 | "y1TErr.="<<errs[1]<< | |
490 | "z0TErr.="<<errs[2]<< | |
491 | "z1TErr.="<<errs[3]<< | |
492 | "y0OErr.="<<errs[4]<< | |
493 | "y1OErr.="<<errs[5]<< | |
494 | "z0OErr.="<<errs[6]<< | |
495 | "z1OErr.="<<errs[7]<< | |
496 | "y0T2Err.="<<errs[8]<< | |
497 | "y1T2Err.="<<errs[9]<< | |
498 | "z0T2Err.="<<errs[10]<< | |
499 | "z1T2Err.="<<errs[11]<< | |
500 | "y0O2Err.="<<errs[12]<< | |
501 | "y1O2Err.="<<errs[13]<< | |
502 | "z0O2Err.="<<errs[14]<< | |
503 | "z1O2Err.="<<errs[15]<< | |
504 | "chi2.="<<&chi2N<< | |
505 | "\n"; | |
506 | } | |
507 | // | |
508 | // delta alpha y0 z0 ky kz | |
509 | Double_t x[6]={0,0,0,0,0,0}; | |
510 | x[1]=fCurrentAlpha; | |
511 | x[2]=(*params[0])[0]; | |
512 | x[3]=(*params[2])[0]; | |
513 | x[4]=(*params[0])[1]; | |
514 | x[5]=(*params[2])[1]; | |
515 | // | |
516 | //Delta y | |
517 | // | |
518 | x[0]= (*params[1])[0]-(*params[0])[0]; | |
519 | fLinearTrackDelta[0]->Fill(x); | |
520 | x[0]/=TMath::Sqrt((*errs[1])[0]*(*errs[1])[0]+(*errs[0])[0]*(*errs[0])[0]); | |
521 | fLinearTrackPull[0]->Fill(x); | |
522 | // | |
523 | //Delta ky | |
524 | // | |
525 | x[0]= 100*((*params[1])[1]-(*params[0])[1]); | |
526 | fLinearTrackDelta[1]->Fill(x); | |
527 | x[0]/=100*TMath::Sqrt((*errs[1])[1]*(*errs[1])[1]+(*errs[0])[1]*(*errs[0])[1]); | |
528 | fLinearTrackPull[1]->Fill(x); | |
529 | // | |
530 | // ky2_0 | |
531 | // | |
532 | x[0]= 100*100*((*params[8])[2]); | |
533 | fLinearTrackDelta[2]->Fill(x); | |
534 | x[0]/=100*100*TMath::Sqrt((*errs[8])[2]*(*errs[8])[2]); | |
535 | fLinearTrackPull[2]->Fill(x); | |
536 | // | |
537 | // ky2_1 | |
538 | // | |
539 | x[0]= 100*100*((*params[9])[2]); | |
540 | fLinearTrackDelta[3]->Fill(x); | |
541 | x[0]/=100*100*TMath::Sqrt((*errs[9])[2]*(*errs[9])[2]); | |
542 | fLinearTrackPull[3]->Fill(x); | |
543 | // | |
544 | // | |
545 | //Delta z | |
546 | // | |
547 | x[0]= (*params[3])[0]-(*params[2])[0]; | |
548 | fLinearTrackDelta[4]->Fill(x); | |
549 | x[0]/=TMath::Sqrt((*errs[3])[0]*(*errs[3])[0]+(*errs[2])[0]*(*errs[2])[0]); | |
550 | fLinearTrackPull[4]->Fill(x); | |
551 | // | |
552 | //Delta kz | |
553 | // | |
554 | x[0]= 100*((*params[3])[1]-(*params[2])[1]); | |
555 | fLinearTrackDelta[5]->Fill(x); | |
556 | x[0]/=100*TMath::Sqrt((*errs[3])[1]*(*errs[3])[1]+(*errs[2])[1]*(*errs[2])[1]); | |
557 | fLinearTrackPull[5]->Fill(x); | |
558 | // | |
559 | // kz2_0 | |
560 | // | |
561 | x[0]= 100*100*((*params[10])[2]); | |
562 | fLinearTrackDelta[6]->Fill(x); | |
563 | x[0]/=100*100*TMath::Sqrt((*errs[10])[2]*(*errs[10])[2]); | |
564 | fLinearTrackPull[6]->Fill(x); | |
565 | // | |
566 | // kz2_1 | |
567 | // | |
568 | x[0]= 100*100*((*params[11])[2]); | |
569 | fLinearTrackDelta[7]->Fill(x); | |
570 | x[0]/=100*100*TMath::Sqrt((*errs[11])[2]*(*errs[11])[2]); | |
571 | fLinearTrackPull[7]->Fill(x); | |
572 | ||
c64c82db | 573 | // |
574 | // rms of track | |
575 | // | |
576 | x[0]= chi2N[0]; | |
577 | fLinearTrackDelta[8]->Fill(x); | |
578 | x[0]= chi2N[1]; | |
579 | fLinearTrackDelta[9]->Fill(x); | |
580 | x[0]= chi2N[2]; | |
581 | fLinearTrackDelta[10]->Fill(x); | |
582 | x[0]= chi2N[3]; | |
583 | fLinearTrackDelta[11]->Fill(x); | |
584 | // | |
0bfb86a0 | 585 | |
586 | ||
587 | for (Int_t ifit=0; ifit<8;ifit++){ | |
588 | delete fitters[ifit]; | |
589 | delete params[ifit]; | |
590 | delete fitters[ifit+8]; | |
591 | delete params[ifit+8]; | |
592 | delete errs[ifit]; | |
593 | delete errs[ifit+8]; | |
594 | } | |
595 | } | |
596 | ||
597 | ||
598 | void AliTPCkalmanFit::Propagate(TTreeSRedirector */*debug*/){ | |
599 | // | |
600 | // Propagate the Kalman | |
601 | // | |
602 | } | |
603 | ||
2497be21 | 604 | void AliTPCkalmanFit::AddCovariance(const char * varName, Double_t sigma){ |
605 | // | |
606 | // | |
607 | // | |
608 | if (fCalibCovar) return; | |
609 | if (!fCalibration) return; | |
610 | if (!fCalibration->FindObject(varName)) return; | |
611 | Int_t ncalibs = fCalibration->GetEntries(); | |
612 | TString strVar(varName); | |
613 | for (Int_t icalib=0;icalib<ncalibs; icalib++){ | |
614 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
615 | if (strVar.CompareTo(transform->GetName())==0){ | |
616 | (*fCalibCovar)(icalib,icalib)+=sigma*sigma; | |
617 | } | |
618 | } | |
619 | } | |
620 | ||
0bfb86a0 | 621 | |
622 | void AliTPCkalmanFit::PropagateTime(Int_t time){ | |
623 | // | |
624 | // Propagate the calibration in time | |
625 | // - Increase covariance matrix | |
626 | // | |
7dfa0700 | 627 | if (!fCalibCovar) return; |
0bfb86a0 | 628 | Int_t ncalibs = fCalibration->GetEntries(); |
629 | Double_t deltaT = 0; | |
630 | if (fLastTimeStamp>0) deltaT = (fLastTimeStamp-time)/(60*60.); // delta T in hours | |
631 | fLastTimeStamp = time; | |
632 | for (Int_t icalib=0;icalib<ncalibs; icalib++){ | |
633 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
98a4cc77 | 634 | if ((*fCalibCovar)(icalib,icalib)<transform->GetSigmaMax()*transform->GetSigmaMax()) |
635 | (*fCalibCovar)(icalib,icalib)+= transform->GetSigma2Time()*TMath::Abs(deltaT); | |
0bfb86a0 | 636 | } |
637 | } | |
638 | ||
639 | ||
640 | void AliTPCkalmanFit::UpdateLinear(AliTrackPoint &point, TTreeSRedirector *debug){ | |
641 | // | |
642 | // Update Kalman | |
643 | // | |
644 | // | |
645 | Int_t ncalibs = fCalibration->GetEntries(); | |
646 | Int_t kNmeas = 2; | |
647 | Int_t nelem = ncalibs+4; | |
648 | TMatrixD &vecXk=*fLinearParam; // X vector | |
649 | TMatrixD &covXk=*fLinearCovar; // X covariance | |
650 | // | |
651 | TMatrixD matHk(kNmeas,nelem); // vector to mesurement | |
652 | TMatrixD vecYk(kNmeas,1); // Innovation or measurement residual | |
653 | TMatrixD vecZk(kNmeas,1); // Innovation or measurement residual | |
654 | TMatrixD measR(kNmeas,kNmeas); | |
655 | TMatrixD matHkT(nelem,kNmeas); // helper matrix Hk transpose | |
656 | TMatrixD matSk(kNmeas,kNmeas); // Innovation (or residual) covariance | |
657 | TMatrixD matKk(nelem,kNmeas); // Optimal Kalman gain | |
658 | TMatrixD covXk2(nelem,nelem); // helper matrix | |
659 | TMatrixD covXk3(nelem,nelem); // helper matrix | |
660 | TMatrixD mat1(nelem,nelem); | |
661 | // | |
662 | // reset matHk | |
663 | for (Int_t iel=0;iel<nelem;iel++){ | |
664 | for (Int_t ip=0;ip<kNmeas;ip++) { | |
665 | matHk(ip,iel)=0; | |
666 | } | |
667 | } | |
668 | for (Int_t iel0=0;iel0<nelem;iel0++) | |
669 | for (Int_t iel1=0;iel1<nelem;iel1++){ | |
670 | if (iel0!=iel1) mat1(iel0,iel1)=0; | |
671 | if (iel0==iel1) mat1(iel0,iel1)=1; | |
672 | } | |
673 | // | |
674 | // fill matrix Hk | |
675 | // | |
676 | Int_t volId = point.GetVolumeID(); | |
677 | Double_t gxyz[3]={point.GetX(), point.GetY(),point.GetZ()}; | |
678 | Double_t rxyz[3]={fCA*gxyz[0]+fSA*gxyz[1],-fSA*gxyz[0]+fCA*gxyz[1] ,point.GetZ()}; | |
679 | // | |
680 | Double_t dxdydz[3]={0,0,0}; | |
681 | Double_t rdxdydz[3]={0,0,0}; | |
682 | ||
683 | for (Int_t icalib=0;icalib<ncalibs;icalib++){ | |
684 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
685 | dxdydz[0] = transform->GetDeltaXYZ(0,volId, 1, gxyz[0], gxyz[1], gxyz[2]); | |
686 | dxdydz[1] = transform->GetDeltaXYZ(1,volId, 1, gxyz[0], gxyz[1], gxyz[2]); | |
687 | dxdydz[2] = transform->GetDeltaXYZ(2,volId, 1, gxyz[0], gxyz[1], gxyz[2]); | |
688 | rdxdydz[0] = fCA*dxdydz[0]+fSA*dxdydz[1]; | |
689 | rdxdydz[1] = -fSA*dxdydz[0]+fCA*dxdydz[1]; | |
690 | rdxdydz[2] = dxdydz[2]; | |
691 | // | |
692 | matHk(0,icalib)= rdxdydz[1]-rdxdydz[0]*(*fLinearParam)(ncalibs+1,0); // shift y + shift x * angleY | |
693 | matHk(1,icalib)= rdxdydz[2]-rdxdydz[0]*(*fLinearParam)(ncalibs+3,0); // shift z + shift x * angleZ | |
694 | } | |
695 | matHk(0,ncalibs+0)=1; | |
696 | matHk(0,ncalibs+1)=rxyz[0]; | |
697 | matHk(1,ncalibs+2)=1; | |
698 | matHk(1,ncalibs+3)=rxyz[0]; | |
699 | // | |
700 | // | |
701 | // | |
702 | vecZk(0,0) = rxyz[1]; | |
703 | vecZk(1,0) = rxyz[2]; | |
704 | measR(0,0) = point.GetCov()[1]; measR(0,1)=0; | |
705 | measR(1,1) = point.GetCov()[2]; measR(1,0)=0; | |
c64c82db | 706 | // |
0bfb86a0 | 707 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual |
708 | matHkT=matHk.T(); matHk.T(); | |
709 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
710 | matSk.Invert(); | |
711 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
712 | // | |
713 | covXk2= (mat1-(matKk*matHk)); | |
714 | covXk3 = covXk2*covXk; | |
715 | // | |
716 | CheckCovariance(covXk3,100); | |
717 | vecXk += matKk*vecYk; // updated vector | |
718 | covXk = covXk3; | |
719 | if (debug){ // dump debug info | |
720 | (*debug)<<"updateLinear"<< | |
721 | // | |
722 | "point.="<<&point<< | |
723 | "vecYk.="<<&vecYk<< | |
724 | "vecZk.="<<&vecZk<< | |
725 | "measR.="<<&measR<< | |
726 | "matHk.="<<&matHk<< | |
727 | "matHkT.="<<&matHkT<< | |
728 | "matSk.="<<&matSk<< | |
729 | "matKk.="<<&matKk<< | |
730 | "covXk2.="<<&covXk2<< | |
731 | "covXk.="<<&covXk<< | |
732 | "vecXk.="<<&vecXk<< | |
733 | "\n"; | |
734 | } | |
735 | } | |
736 | ||
737 | ||
ddf81107 | 738 | AliTrackPointArray * AliTPCkalmanFit::SortPoints(AliTrackPointArray &points){ |
739 | // | |
740 | //Creates the array - points sorted according radius - neccessay for kalman fit | |
741 | // | |
742 | // | |
743 | // 0. choose the frame - rotation angle | |
744 | // | |
745 | Int_t npoints = points.GetNPoints(); | |
746 | if (npoints<1) return 0; | |
747 | Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]); | |
748 | Double_t ca = TMath::Cos(currentAlpha); | |
749 | Double_t sa = TMath::Sin(currentAlpha); | |
750 | // | |
751 | // 1. sort the points | |
752 | // | |
753 | Double_t *rxvector = new Double_t[npoints]; | |
754 | Int_t *indexes = new Int_t[npoints]; | |
755 | for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){ | |
756 | rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint]; | |
757 | } | |
758 | TMath::Sort(npoints, rxvector,indexes,kFALSE); | |
759 | AliTrackPoint point; | |
760 | AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints); | |
761 | for (Int_t ipoint=0; ipoint<npoints; ipoint++){ | |
762 | if (!points.GetPoint(point,indexes[ipoint])) continue; | |
763 | pointsSorted->AddPoint(ipoint,&point); | |
764 | } | |
765 | delete [] rxvector; | |
766 | delete [] indexes; | |
767 | return pointsSorted; | |
768 | } | |
769 | ||
770 | ||
771 | ||
0bfb86a0 | 772 | AliTrackPointArray * AliTPCkalmanFit::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){ |
773 | // | |
774 | // | |
775 | // | |
776 | AliTrackPointArray array(500); | |
777 | Float_t cov[10]; // dummy covariance | |
778 | Int_t npoints=0; | |
779 | for (Int_t i=0;i<6;i++) cov[i]=0.001; | |
780 | for (Int_t i=0;i<500;i++){ | |
30dea41e | 781 | AliTrackPoint point(0, 0, 0, cov, 0,0,0); |
0bfb86a0 | 782 | array.AddPoint(npoints, &point); |
783 | npoints++; | |
784 | } | |
785 | npoints=0; | |
786 | for (Float_t ir = -245; ir<245; ir++){ | |
787 | // | |
788 | // | |
789 | if (TMath::Abs(ir)<80) continue; | |
790 | Double_t ca = TMath::Cos(alpha); | |
791 | Double_t sa = TMath::Sin(alpha); | |
792 | Double_t lx = ir; | |
793 | Double_t ly = y0+lx*ky+gRandom->Gaus(0,err); | |
794 | Double_t lz = z0+lx*kz+gRandom->Gaus(0,err); | |
795 | Double_t gx = lx*ca-ly*sa; | |
796 | Double_t gy = lx*sa+ly*ca; | |
797 | Double_t gz = lz; | |
798 | Double_t galpha = TMath::ATan2(gy,gx); | |
799 | Int_t isec = TMath::Nint((9*galpha/TMath::Pi()+0.5)); | |
800 | if (isec<0) isec+=18; | |
801 | if (gz<0) isec+=18; | |
802 | if (ir>130) isec+=36; | |
803 | // | |
804 | AliTrackPoint point(gx, gy, gz, cov, isec,0,0); | |
805 | array.AddPoint(npoints, &point); | |
806 | npoints++; | |
807 | } | |
808 | AliTrackPointArray *parray = new AliTrackPointArray(npoints); | |
809 | AliTrackPoint point; | |
810 | for (Int_t i=0;i<npoints;i++){ | |
811 | array.GetPoint(point,i); | |
812 | parray->AddPoint(i,&point); | |
813 | } | |
814 | ||
815 | return parray; | |
816 | } | |
817 | ||
818 | void AliTPCkalmanFit::AddCalibration(AliTPCTransformation * calib){ | |
819 | // | |
820 | // Add calibration | |
821 | // | |
822 | if (!fCalibration) fCalibration = new TObjArray(2000); | |
823 | fCalibration->AddLast(calib); | |
824 | } | |
825 | ||
7dfa0700 | 826 | Int_t AliTPCkalmanFit::GetTransformationIndex(const char * trName){ |
827 | // | |
828 | // | |
829 | // | |
830 | if (!fCalibration) return -1; | |
831 | if (!fCalibration->FindObject(trName)) return -1; | |
832 | // | |
833 | Int_t ncalibs = fCalibration->GetEntries(); | |
834 | TString strVar(trName); | |
835 | for (Int_t icalib=0;icalib<ncalibs; icalib++){ | |
836 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
837 | if (strVar.CompareTo(transform->GetName())==0){ | |
838 | return icalib; | |
839 | } | |
840 | } | |
841 | return -1; | |
842 | } | |
843 | ||
844 | ||
845 | ||
0bfb86a0 | 846 | void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){ |
847 | // | |
848 | // | |
849 | // | |
850 | if (!fCalibration) return; | |
851 | Int_t ncalibs = fCalibration->GetEntries(); | |
852 | if (ncalibs==0) return; | |
853 | Int_t npoints = array->GetNPoints(); | |
854 | for (Int_t ipoint=0; ipoint<npoints; ipoint++){ | |
855 | Int_t volId = array->GetVolumeID()[ipoint]; | |
856 | Double_t xyz[3]={array->GetX()[ipoint], array->GetY()[ipoint],array->GetZ()[ipoint]}; | |
857 | Double_t dxdydz[3]={0,0,0}; | |
858 | for (Int_t icalib=0; icalib<ncalibs; icalib++){ | |
859 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
98a4cc77 | 860 | dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]); |
861 | dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]); | |
862 | dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]); | |
0bfb86a0 | 863 | } |
864 | ((Float_t*)array->GetX())[ipoint]+=csign*dxdydz[0]; | |
865 | ((Float_t*)array->GetY())[ipoint]+=csign*dxdydz[1]; | |
866 | ((Float_t*)array->GetZ())[ipoint]+=csign*dxdydz[2]; | |
867 | } | |
868 | } | |
869 | ||
c64c82db | 870 | Bool_t AliTPCkalmanFit::DumpCorelation(Double_t threshold, const char *mask0, const char *mask1){ |
0bfb86a0 | 871 | // |
872 | // | |
873 | // | |
874 | TMatrixD &mat = *fCalibCovar; | |
875 | Int_t nrow= mat.GetNrows(); | |
876 | for (Int_t irow=0; irow<nrow; irow++){ | |
c64c82db | 877 | AliTPCTransformation * trans0 = GetTransformation(irow); |
878 | TString strName0(trans0->GetName()); | |
879 | if (mask0){ | |
880 | if (!strName0.Contains(mask0)) continue; | |
881 | } | |
0bfb86a0 | 882 | for (Int_t icol=irow+1; icol<nrow; icol++){ |
c64c82db | 883 | AliTPCTransformation * trans1 = GetTransformation(icol); |
884 | TString strName1(trans1->GetName()); | |
885 | if (mask1){ | |
886 | if (!strName1.Contains(mask1)) continue; | |
887 | } | |
888 | // | |
0bfb86a0 | 889 | Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol))); |
890 | if (diag<=0){ | |
891 | printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol)); | |
892 | continue; | |
893 | } | |
894 | Double_t corr0 = mat(irow,icol)/diag; | |
895 | if (TMath::Abs(corr0)>threshold){ | |
0bfb86a0 | 896 | printf("%d\t%d\t%s\t%s\t%f\t%f\t%f\n", irow,icol, trans0->GetName(), trans1->GetName(), |
897 | TMath::Sqrt(mat(irow,irow)), TMath::Sqrt(mat(icol,icol)), corr0); | |
898 | } | |
899 | } | |
900 | } | |
2497be21 | 901 | return (nrow>0); |
0bfb86a0 | 902 | } |
903 | ||
2497be21 | 904 | Bool_t AliTPCkalmanFit::DumpCalib(const char *mask, Float_t correlationCut){ |
0bfb86a0 | 905 | // |
906 | // Print calibration entries - name, value, error | |
907 | // | |
908 | TMatrixD &mat = *fCalibCovar; | |
909 | Int_t nrow= mat.GetNrows(); | |
0e9efcbe | 910 | TString strMask(mask); |
2497be21 | 911 | TVectorD vecCorrSum(nrow); |
912 | for (Int_t irow=0; irow<nrow; irow++){ | |
913 | vecCorrSum[irow]=0; | |
914 | for (Int_t icol=0; icol<nrow; icol++){ | |
7dfa0700 | 915 | if (icol==irow) continue; |
2497be21 | 916 | Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol))); |
917 | Double_t corr0 = mat(irow,icol)/diag; | |
918 | vecCorrSum[irow]+=TMath::Abs(corr0); | |
919 | } | |
2497be21 | 920 | vecCorrSum[irow]*=0.5; |
921 | } | |
922 | ||
0bfb86a0 | 923 | for (Int_t irow=0; irow<nrow; irow++){ |
924 | AliTPCTransformation * trans0 = GetTransformation(irow); | |
0e9efcbe | 925 | TString strName(trans0->GetName()); |
926 | if (mask){ | |
927 | if (!strName.Contains(mask)) continue; | |
928 | } | |
2497be21 | 929 | if (vecCorrSum[irow]<correlationCut) continue; |
930 | printf("%d\t%s\t%f\t%f\t%f\n", | |
0bfb86a0 | 931 | irow, |
932 | trans0->GetName(), | |
933 | (*fCalibParam)(irow,0), | |
2497be21 | 934 | TMath::Sqrt(mat(irow,irow)), vecCorrSum[irow]); |
0bfb86a0 | 935 | } |
936 | return (nrow>0); | |
937 | } | |
938 | ||
939 | ||
940 | Bool_t AliTPCkalmanFit::CheckCovariance(TMatrixD &mat, Float_t /*maxEl*/){ | |
941 | // | |
942 | // Check consistency of covariance matrix | |
943 | // + symetrize coavariance matrix | |
944 | Bool_t isOK=kTRUE; | |
945 | Int_t nrow= mat.GetNrows(); | |
946 | for (Int_t irow=0; irow<nrow; irow++){ | |
947 | if (mat(irow,irow)<=0){ | |
948 | printf("Negative covariance\t%d\t%f\n",irow,mat(irow,irow)); | |
949 | isOK=kFALSE; | |
950 | } | |
951 | // if (mat(irow,irow)>maxEl){ | |
952 | // printf("Too big covariance\t%d\t%f\n",irow,mat(irow,irow)); | |
953 | // isOK=kFALSE; | |
954 | // } | |
955 | ||
956 | for (Int_t icol=0; icol<nrow; icol++){ | |
957 | // if (mat(irow,irow) | |
958 | Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol))); | |
959 | if (diag<=0){ | |
960 | printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol)); | |
961 | isOK=kFALSE; | |
962 | continue; | |
963 | } | |
964 | Double_t cov0 = mat(irow,icol)/diag; | |
965 | Double_t cov1 = mat(icol,irow)/diag; | |
966 | if (TMath::Abs(cov0)>1 || TMath::Abs(cov1)>1 ){ | |
967 | printf("Covariance Problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1); | |
968 | isOK=kFALSE; | |
969 | } | |
970 | if (TMath::Abs(cov0-cov1)>0.0000001){ | |
971 | printf("Asymetry problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1); | |
972 | isOK=kFALSE; | |
973 | } | |
974 | // | |
975 | // symetrize the covariance matrix | |
976 | Double_t mean = (mat(irow,icol)+ mat(icol,irow))*0.5; | |
977 | mat(irow,icol)=mean; | |
978 | mat(icol,irow)=mean; | |
979 | } | |
980 | } | |
981 | return isOK; | |
982 | } | |
983 | ||
984 | ||
985 | AliTPCkalmanFit * AliTPCkalmanFit::Test(Int_t ntracks){ | |
986 | // | |
987 | // This is test example | |
988 | // | |
989 | ||
990 | // | |
991 | // 1. Setup transformation | |
992 | // | |
993 | ||
994 | ||
995 | TVectorD fpar(10); | |
996 | AliTPCTransformation * transformation=0; | |
997 | AliTPCkalmanFit * kalmanFit0 = new AliTPCkalmanFit; | |
998 | AliTPCkalmanFit * kalmanFit2 = new AliTPCkalmanFit; | |
999 | // | |
1000 | // Radial scaling | |
1001 | // | |
1002 | for (Int_t iside=0; iside<=1; iside++) | |
1003 | for (Int_t ipar0=0; ipar0<3; ipar0++) | |
1004 | for (Int_t ipar1=0; ipar1<3; ipar1++){ | |
1005 | fpar[0]=ipar0; | |
1006 | fpar[1]=ipar1; | |
1007 | if (ipar0+ipar1==0) continue; | |
1008 | Double_t param = (gRandom->Rndm()-0.5)*0.5; // generate random parameters | |
1009 | char tname[100]; | |
1010 | sprintf(tname,"tscalingR%d%dSide%d",ipar0,ipar1,iside); | |
1011 | transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1); | |
1012 | transformation->SetParams(0,5*0.25,0,&fpar); | |
1013 | kalmanFit0->AddCalibration(transformation); | |
1014 | transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1); | |
1015 | transformation->SetParams(param,0.25,0,&fpar); | |
1016 | kalmanFit2->AddCalibration(transformation); | |
1017 | } | |
1018 | ||
1019 | ||
1020 | // | |
1021 | // 2. Init transformation | |
1022 | // | |
1023 | kalmanFit2->Init(); | |
1024 | kalmanFit0->Init(); | |
1025 | ||
1026 | // | |
1027 | // 3. Run kalman filter | |
1028 | // | |
1029 | Int_t ncalibs = kalmanFit0->fCalibration->GetEntries(); | |
1030 | TVectorD err(ncalibs); | |
1031 | TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTest.root"); | |
1032 | for (Int_t i=0;i<ntracks;i++){ | |
1033 | if (i%100==0) printf("%d\n",i); | |
1034 | Double_t alpha = gRandom->Rndm()*TMath::TwoPi(); | |
1035 | Double_t y0 = (gRandom->Rndm()-0.5)*180; | |
1036 | Double_t z0 = (gRandom->Rndm()-0.5)*250*2; | |
1037 | Double_t ky = (gRandom->Rndm()-0.5)*1; | |
1038 | Double_t kz = (gRandom->Rndm()-0.5)*1; | |
1039 | //generate random TPC track | |
1040 | AliTrackPointArray * array = AliTPCkalmanFit::MakePointArrayLinear(alpha,y0,z0, ky, kz,0.04); | |
1041 | AliTrackPointArray * arrayB = new AliTrackPointArray(*array); // backup | |
1042 | kalmanFit2->ApplyCalibration(array,1.); // misalign ideal track | |
1043 | for (Int_t icalib=0; icalib<ncalibs; icalib++){ | |
1044 | err[icalib] = TMath::Sqrt((*kalmanFit0->fCalibCovar)(icalib,icalib)); | |
1045 | } | |
1046 | // | |
1047 | (*pcstream)<<"dump0"<< | |
1048 | "alpha="<<alpha<< | |
1049 | "y0="<<y0<< | |
1050 | "z0="<<z0<< | |
1051 | "ky="<<ky<< | |
1052 | "lz="<<kz<< | |
1053 | "p.="<<array<< | |
1054 | "pB.="<<arrayB<< | |
1055 | "cparam.="<<kalmanFit0->fCalibParam<< | |
1056 | "ccovar.="<<kalmanFit0->fCalibCovar<< | |
1057 | "err.="<<&err<< | |
1058 | "gparam.="<<kalmanFit2->fCalibParam<< | |
1059 | "gcovar.="<<kalmanFit2->fCalibCovar<< | |
1060 | "\n"; | |
1061 | if (i%20==0) { | |
c64c82db | 1062 | kalmanFit0->FitTrackLinear(*array,pcstream); // fit track - dump intermediate results |
0bfb86a0 | 1063 | }else{ |
c64c82db | 1064 | kalmanFit0->FitTrackLinear(*array,0); // fit track + calibration |
0bfb86a0 | 1065 | } |
1066 | kalmanFit0->DumpTrackLinear(*array,pcstream); // dump track residuals to the tree + fill histograms | |
1067 | } | |
1068 | pcstream->GetFile()->cd(); | |
1069 | kalmanFit0->Write("kalmanFit"); | |
1070 | delete pcstream; | |
1071 | return kalmanFit0; | |
1072 | } | |
1073 | ||
1074 | ||
7dfa0700 | 1075 | // Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){ |
1076 | // // | |
1077 | // // function for visualization purposes | |
1078 | // // | |
1079 | // if (!fCalibration) return 0; | |
1080 | // Int_t ncalibs = fCalibration->GetEntries(); | |
1081 | // if (ncalibs==0) return 0; | |
1082 | // Double_t dxdydz[3]={0,0,0}; | |
1083 | // // | |
1084 | // if (volID<0){ | |
1085 | // Double_t alpha = TMath::ATan2(y,x); | |
1086 | // Double_t r = TMath::Sqrt(y*y+x*x); | |
1087 | // volID = TMath::Nint(9*alpha/TMath::Pi()-0.5); | |
1088 | // if (volID<0) volID+=18; | |
1089 | // if (z<0) volID+=18; | |
1090 | // if (r>120) volID+=36; | |
1091 | // } | |
1092 | // for (Int_t icalib=0; icalib<ncalibs; icalib++){ | |
1093 | // AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
1094 | // Double_t param = (*fCalibParam)(icalib,0); | |
1095 | // dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z); | |
1096 | // } | |
1097 | // return dxdydz[coord]; | |
1098 | // } | |
1099 | ||
1100 | // Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){ | |
1101 | // // | |
1102 | // // | |
1103 | // // | |
1104 | // if (AliTPCkalmanFit::fgInstance==0) return 0; | |
1105 | // return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z); | |
1106 | // } | |
1107 | ||
1108 | ||
1109 | ||
1110 | ||
1111 | ||
1112 | Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){ | |
0bfb86a0 | 1113 | // |
1114 | // function for visualization purposes | |
1115 | // | |
7dfa0700 | 1116 | // coord - coordinate for output |
1117 | // - 0 -X | |
1118 | // 1 -Y | |
1119 | // 2 -Z | |
1120 | // 3 -R | |
1121 | // 4 -RPhi | |
1122 | // 5 -Z | |
1123 | // | |
1124 | //icoordsys - type of coordinate system for input | |
1125 | // - 0 - x,y,z | |
1126 | // - 1 - r,phi,z | |
1127 | // | |
0bfb86a0 | 1128 | if (!fCalibration) return 0; |
1129 | Int_t ncalibs = fCalibration->GetEntries(); | |
1130 | if (ncalibs==0) return 0; | |
7dfa0700 | 1131 | Double_t xyz[3]={0,0,0}; |
1132 | Double_t dxdydz[6]={0,0,0,0,0,0}; | |
1b554eb4 | 1133 | Double_t alpha=0; |
1134 | Double_t r=0; | |
f9f69968 | 1135 | if(icoordsys==0){alpha=TMath::ATan2(y,x); r =TMath::Sqrt(y*y+x*x);} |
1136 | if(icoordsys==1){alpha=y; r=x;} | |
7dfa0700 | 1137 | Double_t ca = TMath::Cos(alpha); |
1138 | Double_t sa = TMath::Sin(alpha); | |
f9f69968 | 1139 | if(icoordsys==0){xyz[0]=x; xyz[1]=y; xyz[2]=z;} |
1140 | if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;} | |
0bfb86a0 | 1141 | // |
1142 | if (volID<0){ | |
7dfa0700 | 1143 | // Double_t alpha = TMath::ATan2(y,x); |
1144 | //Double_t r = TMath::Sqrt(y*y+x*x); | |
0bfb86a0 | 1145 | volID = TMath::Nint(9*alpha/TMath::Pi()-0.5); |
1146 | if (volID<0) volID+=18; | |
1147 | if (z<0) volID+=18; | |
1148 | if (r>120) volID+=36; | |
1149 | } | |
1150 | for (Int_t icalib=0; icalib<ncalibs; icalib++){ | |
1151 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
1152 | Double_t param = (*fCalibParam)(icalib,0); | |
7dfa0700 | 1153 | for (Int_t icoord=0;icoord<6;icoord++){ |
1154 | dxdydz[icoord] += transform->GetDeltaXYZ(icoord,volID, param, xyz[0],xyz[1],xyz[2]); | |
1155 | } | |
0bfb86a0 | 1156 | } |
7dfa0700 | 1157 | |
0bfb86a0 | 1158 | return dxdydz[coord]; |
1159 | } | |
1160 | ||
7dfa0700 | 1161 | |
1162 | Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){ | |
0bfb86a0 | 1163 | // |
1164 | // | |
1165 | // | |
1166 | if (AliTPCkalmanFit::fgInstance==0) return 0; | |
7dfa0700 | 1167 | return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID, icoordsys,x,y,z); |
0bfb86a0 | 1168 | } |
1169 | ||
7dfa0700 | 1170 | |
1171 | ||
1172 | ||
1173 | ||
1174 | Double_t AliTPCkalmanFit::GetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){ | |
1175 | ||
1176 | Int_t ncalibs = fCalibration->GetEntries(); | |
e429bd5a | 1177 | if (calibID>=ncalibs) return 0; |
7dfa0700 | 1178 | //Int_t volID=-1; |
1179 | //Double_t xyz[3]={x,y,z}; | |
1b554eb4 | 1180 | Double_t r=0; |
1181 | Double_t alpha=0; | |
e0d8af3e | 1182 | if(icoordsys==0){r=TMath::Sqrt(x*x+y*y); alpha = TMath::ATan2(y,x);} |
1183 | if(icoordsys==1){r=x; alpha = y;} | |
7dfa0700 | 1184 | Double_t ca = TMath::Cos(alpha); |
1185 | Double_t sa = TMath::Sin(alpha); | |
1b554eb4 | 1186 | Double_t xyz[3]={0,0,0}; |
f9f69968 | 1187 | if(icoordsys==0){xyz[0]=x;xyz[1]=y;xyz[2]=z;} |
1188 | if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;} | |
7dfa0700 | 1189 | //xyz[3]=param; xyz[4]=volID; |
1190 | ||
1191 | if (volID<0){ | |
1192 | //Double_t alpha = TMath::ATan2(xyz[1],xyz[0]); | |
1193 | //Double_t r = TMath::Sqrt(xyz[1]*xyz[1]+xyz[0]*xyz[0]); | |
1194 | volID = TMath::Nint(9*alpha/TMath::Pi()-0.5); | |
1195 | if (volID<0) volID+=18; | |
1196 | if (xyz[2]<0) volID+=18; | |
1197 | if (r>120) volID+=36; | |
1198 | } | |
1199 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(calibID); | |
1200 | //transform->SetInstance(transform); | |
1201 | Double_t param = (*fCalibParam)(calibID,0); | |
1202 | Double_t delta = (Double_t)transform->GetDeltaXYZ(coord,volID, param, xyz[0],xyz[1],xyz[2]); | |
1203 | ||
1204 | return delta; | |
1205 | } | |
1206 | ||
1207 | Double_t AliTPCkalmanFit::SGetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){ | |
1208 | // | |
1209 | // | |
1210 | // | |
1211 | if (AliTPCkalmanFit::fgInstance==0) return 0; | |
1212 | return AliTPCkalmanFit::fgInstance->GetTPCtransXYZ(coord, volID, calibID,icoordsys,x,y,z); | |
1213 | } | |
e0d8af3e | 1214 | |
1215 | ||
f9f69968 | 1216 | void AliTPCkalmanFit::MakeTreeTrans(TTreeSRedirector *debug, const char *treeName){ |
e0d8af3e | 1217 | // |
1218 | // Make the Tree before and after current calibration | |
1219 | // | |
1220 | if (!fCalibParam) { | |
1221 | AliError("Kalman Fit not initialized"); | |
1222 | return; | |
1223 | } | |
1224 | // | |
1225 | // | |
1226 | // | |
1227 | const Int_t ncalibs = fCalibration->GetEntries(); | |
1228 | TMatrixD dxdydz(ncalibs,5); | |
f9f69968 | 1229 | Double_t * adx = new Double_t[ncalibs]; |
1230 | Double_t * ady = new Double_t[ncalibs]; | |
1231 | Double_t * adz = new Double_t[ncalibs]; | |
1232 | Double_t * adr = new Double_t[ncalibs]; | |
1233 | Double_t * adrphi = new Double_t[ncalibs]; | |
e0d8af3e | 1234 | |
1235 | Double_t x[3]; | |
1236 | for (x[0]=-250.;x[0]<=250.;x[0]+=10.){ | |
1237 | for (x[1]=-250.;x[1]<=250.;x[1]+=10.){ | |
1238 | for (x[2]=-250.;x[2]<=250.;x[2]+=20.) { | |
1239 | Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]); | |
1240 | if (r<20) continue; | |
1241 | if (r>260) continue; | |
1242 | //Double_t z = x[2]; | |
1243 | Double_t phi=TMath::ATan2(x[1],x[0]); | |
1244 | Double_t ca=TMath::Cos(phi); | |
1245 | Double_t sa=TMath::Sin(phi); | |
1246 | Double_t dx=0; | |
1247 | Double_t dy=0; | |
1248 | Double_t dz=0; | |
1249 | Double_t dr=0; | |
1250 | Double_t rdphi=0; | |
1251 | ||
f9f69968 | 1252 | Int_t volID= TMath::Nint(9*phi/TMath::Pi()-0.5); |
1253 | if (volID<0) volID+=18; | |
1254 | if (x[2]<0) volID+=18; //C-side | |
1255 | if (r>120) volID+=36; //outer | |
1256 | Double_t volalpha=(volID+0.5)*TMath::Pi()/9; | |
1257 | Double_t cva=TMath::Cos(volalpha); | |
1258 | Double_t sva=TMath::Sin(volalpha); | |
1259 | ||
1260 | Double_t lx=x[0]*cva+x[1]*sva; | |
1261 | Double_t ly=-x[0]*sva+x[1]*cva; | |
1262 | ||
e0d8af3e | 1263 | |
1264 | for(Int_t icalib=0;icalib<ncalibs;icalib++){ | |
1265 | for(Int_t icoord=0;icoord<5;icoord++){ | |
1266 | dxdydz(icalib,icoord)= GetTPCtransXYZ(icoord, -1, icalib, 0, x[0], x[1], x[2]); | |
1267 | } | |
1268 | } | |
1269 | dx=GetTPCDeltaXYZ(0, -1, 0, x[0], x[1], x[2]); | |
1270 | dy=GetTPCDeltaXYZ(1, -1, 0, x[0], x[1], x[2]); | |
1271 | dz=GetTPCDeltaXYZ(2, -1, 0, x[0], x[1], x[2]); | |
1272 | dr=GetTPCDeltaXYZ(3, -1, 0, x[0], x[1], x[2]); | |
1273 | rdphi=GetTPCDeltaXYZ(4, -1, 0, x[0], x[1], x[2]); | |
1274 | ||
1275 | ||
1276 | if(debug){ | |
f9f69968 | 1277 | TTreeStream &cstream= |
1278 | (*debug)<<treeName<< | |
e0d8af3e | 1279 | "x="<<x[0]<< |
1280 | "y="<<x[1]<< | |
1281 | "z="<<x[2]<< | |
1282 | "r="<<r<< | |
1283 | "ca="<<ca<< | |
1284 | "sa="<<sa<< | |
f9f69968 | 1285 | "lx="<<lx<< |
1286 | "ly="<<ly<< | |
1287 | "sector="<<volID<< | |
1288 | "phi="<<phi<< | |
e0d8af3e | 1289 | "dx="<<dx<< |
1290 | "dy="<<dy<< | |
1291 | "dz="<<dz<< | |
1292 | "dr="<<dr<< | |
1293 | "rdphi="<<rdphi<< | |
f9f69968 | 1294 | "dxdydz.="<<&dxdydz; |
1295 | for (Int_t icalib=0;icalib<ncalibs;icalib++){ | |
1296 | AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib); | |
1297 | char tname[1000]; | |
1298 | // | |
1299 | sprintf(tname,"dx%s=",transform->GetName()); | |
1300 | adx[icalib] =dxdydz(icalib,0); | |
1301 | cstream<<tname<<adx[icalib]; | |
1302 | sprintf(tname,"dy%s=",transform->GetName()); | |
1303 | ady[icalib] =dxdydz(icalib,1); | |
1304 | cstream<<tname<<ady[icalib]; | |
1305 | sprintf(tname,"dz%s=",transform->GetName()); | |
1306 | adz[icalib] =dxdydz(icalib,2); | |
1307 | cstream<<tname<<adz[icalib]; | |
1308 | // | |
1309 | sprintf(tname,"dr%s=",transform->GetName()); | |
1310 | adr[icalib] =dxdydz(icalib,3); | |
1311 | cstream<<tname<<adr[icalib]; | |
1312 | sprintf(tname,"rdphi%s=",transform->GetName()); | |
1313 | adrphi[icalib] =dxdydz(icalib,4); | |
1314 | cstream<<tname<<adrphi[icalib]; | |
1315 | } | |
1316 | cstream<<"\n"; | |
e0d8af3e | 1317 | } |
1318 | } | |
1319 | } | |
1320 | Printf("x0=%f finished",x[0]); | |
1321 | } | |
f9f69968 | 1322 | |
e0d8af3e | 1323 | } |