]>
Commit | Line | Data |
---|---|---|
69209e2c | 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 | // This class stores a tracklet (a track that lives only in a single TPC | |
18 | // sector). Its objects can be constructed out of TPCseeds, that are | |
19 | // holding the necessary cluster information. | |
20 | //// | |
21 | //// | |
22 | //// | |
23 | ||
69209e2c | 24 | #include "AliTPCTracklet.h" |
25 | #include "TObjArray.h" | |
9318a5b4 | 26 | #include "TLinearFitter.h" |
69209e2c | 27 | #include "AliTPCseed.h" |
28 | #include "AliESDVertex.h" | |
9318a5b4 | 29 | #include "AliTracker.h" |
30 | #include "TTreeStream.h" | |
31 | #include "TRandom3.h" | |
32 | #include "TDecompChol.h" | |
33 | ||
34 | #include <iostream> | |
35 | using namespace std; | |
69209e2c | 36 | |
37 | ClassImp(AliTPCTracklet) | |
38 | ||
9318a5b4 | 39 | const Double_t AliTPCTracklet::kB2C=0.299792458e-3; |
967eae0d | 40 | Float_t AliTPCTracklet::fgEdgeCutY=3; |
41 | Float_t AliTPCTracklet::fgEdgeCutX=0; | |
9318a5b4 | 42 | |
69209e2c | 43 | AliTPCTracklet::AliTPCTracklet() |
9318a5b4 | 44 | : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0), |
45 | fInner(0),fPrimary(0) { | |
69209e2c | 46 | //// |
47 | // The default constructor. It is intended to be used for I/O only. | |
48 | //// | |
49 | } | |
50 | ||
9318a5b4 | 51 | AliTPCTracklet::AliTPCTracklet(const AliTPCseed *track,Int_t sector, |
52 | TrackType type,Bool_t storeClusters) | |
53 | : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0), | |
54 | fInner(0),fPrimary(0) { | |
69209e2c | 55 | //// |
56 | // Contructor for a tracklet out of a track. Only clusters within a given | |
57 | // sector are used. | |
58 | /// | |
69209e2c | 59 | |
9318a5b4 | 60 | //TODO: only kalman works |
61 | ||
69209e2c | 62 | for (Int_t i=0;i<160;++i) { |
9318a5b4 | 63 | AliTPCclusterMI *c=track->GetClusterPointer(i); |
967eae0d | 64 | if (c && RejectCluster(c)) continue; |
9318a5b4 | 65 | if (c&&c->GetDetector()==sector) |
66 | ++fNClusters; | |
69209e2c | 67 | } |
9318a5b4 | 68 | |
69 | if (storeClusters) { | |
70 | fClusters=new AliTPCclusterMI[fNClusters]; | |
71 | for (Int_t i=0;i<160;++i) { | |
72 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
967eae0d | 73 | if (c && RejectCluster(c)) continue; |
9318a5b4 | 74 | if (c&&c->GetDetector()==sector) |
d47d1dfd | 75 | fClusters[fNStoredClusters]=*c; |
9318a5b4 | 76 | ++fNStoredClusters; |
69209e2c | 77 | } |
78 | } | |
9318a5b4 | 79 | |
80 | switch (type) { | |
81 | case kKalman: | |
82 | FitKalman(track,sector); | |
83 | break; | |
84 | case kLinear: | |
85 | case kQuadratic: | |
86 | FitLinear(track,sector,type); | |
87 | break; | |
88 | case kRiemann: | |
89 | FitRiemann(track,sector); | |
90 | break; | |
69209e2c | 91 | } |
69209e2c | 92 | |
9318a5b4 | 93 | } |
69209e2c | 94 | |
c94a79e1 | 95 | AliTPCTracklet::AliTPCTracklet(const TObjArray &/*clusters*/,Int_t sector, |
96 | TrackType /*type*/,Bool_t /*storeClusters*/) | |
97 | : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0), | |
98 | fInner(0),fPrimary(0) { | |
9318a5b4 | 99 | //TODO: write it! |
69209e2c | 100 | } |
101 | ||
102 | AliTPCTracklet::AliTPCTracklet(const AliTPCTracklet &t) | |
c94a79e1 | 103 | : TObject(t),fNClusters(t.fNClusters),fNStoredClusters(t.fNStoredClusters),fClusters(0), |
9318a5b4 | 104 | fSector(t.fSector),fOuter(0),fInner(0), |
69209e2c | 105 | fPrimary(0) { |
106 | //// | |
107 | // The copy constructor. You can copy tracklets! | |
108 | //// | |
109 | ||
9318a5b4 | 110 | if (t.fClusters) { |
111 | fClusters=new AliTPCclusterMI[t.fNStoredClusters]; | |
112 | for (int i=0;i<t.fNStoredClusters;++i) | |
113 | fClusters[i]=t.fClusters[i]; | |
114 | } | |
69209e2c | 115 | if (t.fOuter) |
116 | fOuter=new AliExternalTrackParam(*t.fOuter); | |
117 | if (t.fInner) | |
118 | fInner=new AliExternalTrackParam(*t.fInner); | |
119 | if (t.fPrimary) | |
120 | fPrimary=new AliExternalTrackParam(*t.fPrimary); | |
121 | } | |
122 | ||
123 | AliTPCTracklet& AliTPCTracklet::operator=(const AliTPCTracklet &t) { | |
124 | //// | |
125 | // The assignment constructor. You can assign tracklets! | |
126 | //// | |
69209e2c | 127 | if (this!=&t) { |
9318a5b4 | 128 | fNClusters=t.fNClusters; |
129 | fNStoredClusters=fNStoredClusters; | |
130 | delete fClusters; | |
131 | if (t.fClusters) { | |
132 | fClusters=new AliTPCclusterMI[t.fNStoredClusters]; | |
133 | for (int i=0;i<t.fNStoredClusters;++i) | |
134 | fClusters[i]=t.fClusters[i]; | |
135 | } | |
136 | else | |
137 | fClusters=0; | |
138 | fSector=t.fSector; | |
69209e2c | 139 | if (t.fOuter) { |
140 | if (fOuter) | |
141 | *fOuter=*t.fOuter; | |
142 | else | |
143 | fOuter=new AliExternalTrackParam(*t.fOuter); | |
144 | } | |
145 | else { | |
146 | delete fOuter; | |
147 | fOuter=0; | |
148 | } | |
149 | ||
150 | if (t.fInner) { | |
151 | if (fInner) | |
152 | *fInner=*t.fInner; | |
153 | else | |
154 | fInner=new AliExternalTrackParam(*t.fInner); | |
155 | } | |
156 | else { | |
157 | delete fInner; | |
158 | fInner=0; | |
159 | } | |
160 | ||
161 | if (t.fPrimary) { | |
162 | if (fPrimary) | |
163 | *fPrimary=*t.fPrimary; | |
164 | else | |
165 | fPrimary=new AliExternalTrackParam(*t.fPrimary); | |
166 | } | |
167 | else { | |
168 | delete fPrimary; | |
169 | fPrimary=0; | |
170 | } | |
171 | } | |
172 | return *this; | |
173 | } | |
174 | ||
175 | AliTPCTracklet::~AliTPCTracklet() { | |
176 | // | |
177 | // The destructor. Yes, you can even destruct tracklets. | |
178 | // | |
9318a5b4 | 179 | delete fClusters; |
69209e2c | 180 | delete fOuter; |
181 | delete fInner; | |
182 | delete fPrimary; | |
183 | } | |
184 | ||
ef0a7e89 | 185 | |
186 | ||
187 | /* | |
9318a5b4 | 188 | void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) { |
8b3c60d8 | 189 | // |
190 | // Fit using Kalman filter | |
191 | // | |
9318a5b4 | 192 | AliTPCseed *t=new AliTPCseed(*track); |
193 | if (!t->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-t->GetAlpha())) { | |
194 | delete t; | |
195 | return; | |
196 | } | |
197 | // fit from inner to outer row | |
ef0a7e89 | 198 | Double_t covar[15]; |
199 | for (Int_t i=0;i<15;i++) covar[i]=0; | |
200 | covar[0]=10.*10.; | |
201 | covar[2]=10.*10.; | |
202 | covar[5]=10.*10./(64.*64.); | |
203 | covar[9]=10.*10./(64.*64.); | |
204 | covar[14]=1*1; | |
205 | ||
206 | // | |
9318a5b4 | 207 | AliTPCseed *outerSeed=new AliTPCseed(*t); |
208 | Int_t n=0; | |
209 | for (Int_t i=0;i<160;++i) { | |
ef0a7e89 | 210 | |
9318a5b4 | 211 | AliTPCclusterMI *c=t->GetClusterPointer(i); |
967eae0d | 212 | if (c && RejectCluster(c,outerSeed)) continue; |
9318a5b4 | 213 | if (c&&c->GetDetector()==sector) { |
214 | if (n==1) { | |
215 | outerSeed->ResetCovariance(100.); | |
ef0a7e89 | 216 | outerSeed->AddCovariance(covar); |
9318a5b4 | 217 | } |
218 | ++n; | |
219 | Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; | |
ef0a7e89 | 220 | Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation |
9318a5b4 | 221 | if (!outerSeed->PropagateTo(r[0]) || |
222 | !static_cast<AliExternalTrackParam*>(outerSeed)->Update(&r[1],cov)) { | |
223 | delete outerSeed; | |
224 | outerSeed=0; | |
225 | break; | |
226 | } | |
227 | } | |
228 | } | |
229 | if (outerSeed) | |
230 | fOuter=new AliExternalTrackParam(*outerSeed); | |
9318a5b4 | 231 | // fit from outer to inner rows |
ef0a7e89 | 232 | // AliTPCseed *innerSeed=new AliTPCseed(*t); |
233 | AliTPCseed *innerSeed=0; | |
234 | if (fOuter) innerSeed=new AliTPCseed(*outerSeed); | |
235 | if (!innerSeed) innerSeed=new AliTPCseed(*t); | |
236 | delete outerSeed; | |
237 | ||
9318a5b4 | 238 | n=0; |
239 | for (Int_t i=159;i>=0;--i) { | |
240 | AliTPCclusterMI *c=t->GetClusterPointer(i); | |
967eae0d | 241 | if (c && RejectCluster(c, innerSeed)) continue; |
9318a5b4 | 242 | if (c&&c->GetDetector()==sector) { |
243 | if (n==1) { | |
244 | innerSeed->ResetCovariance(100.); | |
ef0a7e89 | 245 | innerSeed->AddCovariance(covar); |
9318a5b4 | 246 | } |
247 | ++n; | |
248 | Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; | |
ef0a7e89 | 249 | Double_t cov[3]={0.01,0.,0.01}; |
9318a5b4 | 250 | if (!innerSeed->PropagateTo(r[0]) || |
251 | !static_cast<AliExternalTrackParam*>(innerSeed)->Update(&r[1],cov)) { | |
252 | delete innerSeed; | |
253 | innerSeed=0; | |
254 | break; | |
255 | } | |
256 | } | |
257 | } | |
258 | if (innerSeed) | |
259 | fInner=new AliExternalTrackParam(*innerSeed); | |
260 | // propagate to the primary vertex | |
261 | if (innerSeed) { | |
262 | AliTPCseed *primarySeed=new AliTPCseed(*innerSeed); | |
263 | Double_t pos[]={0.,0.,0.}; | |
264 | Double_t sigma[]={.1,.1,.1}; //TODO: is this correct? | |
265 | AliESDVertex vertex(pos,sigma); | |
266 | if (primarySeed->PropagateToVertex(&vertex)) | |
267 | fPrimary=new AliExternalTrackParam(*primarySeed); | |
268 | delete primarySeed; | |
269 | // for better comparison one does not want to have alpha changed... | |
7eaa723e | 270 | if (fPrimary) if (!fPrimary->Rotate(fInner->GetAlpha())) { |
9318a5b4 | 271 | delete fPrimary; |
272 | fPrimary=0; | |
273 | } | |
274 | } | |
275 | delete innerSeed; | |
276 | ||
277 | delete t; | |
278 | } | |
ef0a7e89 | 279 | */ |
280 | ||
281 | ||
282 | void AliTPCTracklet::FitKalman(const AliTPCseed *seed,Int_t sector) { | |
283 | // | |
284 | // Fit using Kalman filter | |
285 | // | |
286 | AliTPCseed *track=new AliTPCseed(*seed); | |
287 | if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) { | |
288 | delete track; | |
289 | return; | |
290 | } | |
291 | // fit from inner to outer row | |
292 | Double_t covar[15]; | |
293 | for (Int_t i=0;i<15;i++) covar[i]=0; | |
294 | covar[0]=10.*10.; | |
295 | covar[2]=10.*10.; | |
296 | covar[5]=10.*10./(64.*64.); | |
297 | covar[9]=10.*10./(64.*64.); | |
298 | covar[14]=1*1; | |
299 | Float_t xmin=1000, xmax=-10000; | |
300 | Int_t imin=158, imax=0; | |
301 | for (Int_t i=0;i<160;i++) { | |
302 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
303 | if (!c) continue; | |
304 | if (c->GetDetector()!=sector) continue; | |
305 | if (c->GetX()<xmin) xmin=c->GetX(); | |
306 | if (c->GetX()>xmax) xmax=c->GetX(); | |
307 | if (i<imin) imin=i; | |
308 | if (i>imax) imax=i; | |
309 | } | |
310 | if(imax-imin<10) { | |
311 | delete track; | |
312 | return; | |
313 | } | |
314 | ||
315 | for (Float_t x=track->GetX(); x<xmin; x++) track->PropagateTo(x); | |
316 | track->AddCovariance(covar); | |
317 | // | |
318 | AliExternalTrackParam paramIn; | |
319 | AliExternalTrackParam paramOut; | |
320 | Bool_t isOK=kTRUE; | |
321 | // | |
322 | // | |
323 | // | |
324 | for (Int_t i=imin; i<=imax; i++){ | |
325 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
326 | if (!c) continue; | |
ef0a7e89 | 327 | Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; |
328 | Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation | |
91350f5b | 329 | AliTPCseed::GetError(c, track,cov[0],cov[2]); |
330 | cov[0]*=cov[0]; | |
331 | cov[2]*=cov[2]; | |
ef0a7e89 | 332 | if (!track->PropagateTo(r[0])) { |
333 | isOK=kFALSE; | |
334 | break; | |
335 | } | |
42b40d07 | 336 | if (RejectCluster(c,track)) continue; |
ef0a7e89 | 337 | if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE; |
338 | } | |
339 | if (!isOK) { delete track; return;} | |
340 | track->AddCovariance(covar); | |
341 | // | |
342 | // | |
343 | // | |
344 | for (Int_t i=imax; i>=imin; i--){ | |
345 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
346 | if (!c) continue; | |
ef0a7e89 | 347 | Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; |
91350f5b | 348 | Double_t cov[3]={0.01,0.,0.01}; |
349 | AliTPCseed::GetError(c, track,cov[0],cov[2]); | |
350 | cov[0]*=cov[0]; | |
351 | cov[2]*=cov[2]; | |
ef0a7e89 | 352 | if (!track->PropagateTo(r[0])) { |
353 | isOK=kFALSE; | |
354 | break; | |
355 | } | |
42b40d07 | 356 | if (RejectCluster(c,track)) continue; |
ef0a7e89 | 357 | if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE; |
358 | } | |
359 | if (!isOK) { delete track; return;} | |
360 | paramIn = *track; | |
361 | track->AddCovariance(covar); | |
362 | // | |
363 | // | |
364 | for (Int_t i=imin; i<=imax; i++){ | |
365 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
366 | if (!c) continue; | |
ef0a7e89 | 367 | Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; |
91350f5b | 368 | Double_t cov[3]={0.01,0.,0.01}; |
369 | AliTPCseed::GetError(c, track,cov[0],cov[2]); | |
370 | cov[0]*=cov[0]; | |
371 | cov[2]*=cov[2]; | |
ef0a7e89 | 372 | if (!track->PropagateTo(r[0])) { |
373 | isOK=kFALSE; | |
374 | break; | |
375 | } | |
42b40d07 | 376 | if (RejectCluster(c,track)) continue; |
ef0a7e89 | 377 | if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE; |
378 | } | |
379 | if (!isOK) { delete track; return;} | |
380 | paramOut=*track; | |
381 | // | |
382 | // | |
383 | // | |
384 | fOuter=new AliExternalTrackParam(paramOut); | |
385 | fInner=new AliExternalTrackParam(paramIn); | |
386 | // | |
387 | ||
388 | ||
389 | // // propagate to the primary vertex | |
390 | // if (fInner) { | |
391 | // AliExternalTrackParam *param= new AliExternalTrackParam(*fInner); | |
392 | // Double_t pos[]={0.,0.,0.}; | |
393 | // Double_t sigma[]={.1,.1,.1}; //TODO: is this correct? | |
394 | // AliESDVertex vertex(pos,sigma); | |
395 | // if (param->PropagateToVertex(&vertex)) | |
396 | // fPrimary=new AliExternalTrackParam(*param); | |
397 | // delete param; | |
398 | // // for better comparison one does not want to have alpha changed... | |
399 | // if (fPrimary) if (!fPrimary->Rotate(fInner->GetAlpha())) { | |
400 | // delete fPrimary; | |
401 | // fPrimary=0; | |
402 | // } | |
403 | // } | |
404 | ||
405 | delete track; | |
406 | } | |
407 | ||
408 | ||
409 | ||
9318a5b4 | 410 | |
411 | void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector, | |
412 | TrackType type) { | |
413 | TLinearFitter fy(1); | |
414 | TLinearFitter fz(1); | |
415 | fy.StoreData(kFALSE); | |
416 | fz.StoreData(kFALSE); | |
417 | switch (type) { | |
418 | case kLinear: | |
419 | fy.SetFormula("1 ++ x"); | |
420 | fz.SetFormula("1 ++ x"); | |
421 | break; | |
422 | case kQuadratic: | |
423 | fy.SetFormula("1 ++ x ++ x*x"); | |
424 | fz.SetFormula("1 ++ x"); | |
425 | break; | |
c94a79e1 | 426 | case kKalman: |
427 | case kRiemann: | |
428 | break; | |
9318a5b4 | 429 | } |
430 | Double_t xmax=-1.; | |
431 | Double_t xmin=1000.; | |
432 | for (Int_t i=0;i<160;++i) { | |
433 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
967eae0d | 434 | if (c && RejectCluster(c)) continue; |
9318a5b4 | 435 | if (c&&c->GetDetector()==sector) { |
436 | Double_t x=c->GetX(); | |
437 | fy.AddPoint(&x,c->GetY()); | |
438 | fz.AddPoint(&x,c->GetZ()); | |
439 | xmax=TMath::Max(xmax,x); | |
440 | xmin=TMath::Min(xmin,x); | |
441 | } | |
442 | } | |
443 | fy.Eval(); | |
444 | fz.Eval(); | |
445 | Double_t a[3]={fy.GetParameter(0), | |
446 | fy.GetParameter(1), | |
447 | type==kQuadratic?fy.GetParameter(2):0.}; | |
448 | Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0), | |
449 | fy.GetCovarianceMatrixElement(1,0), | |
450 | fy.GetCovarianceMatrixElement(1,1), | |
451 | type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0., | |
452 | type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0., | |
453 | type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.}; | |
454 | for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters; | |
455 | Double_t b[2]={fz.GetParameter(0), | |
456 | fz.GetParameter(1)}; | |
457 | Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0), | |
458 | fz.GetCovarianceMatrixElement(1,0), | |
459 | fz.GetCovarianceMatrixElement(1,1)}; | |
460 | for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters; | |
461 | Double_t p[5]; | |
462 | Double_t c[15]; | |
463 | Double_t alpha=track->GetAlpha(); | |
464 | Quadratic2Helix(a,ca,b,cb,0.,p,c); | |
465 | fPrimary=new AliExternalTrackParam(0.,alpha,p,c); | |
466 | Quadratic2Helix(a,ca,b,cb,xmin,p,c); | |
467 | fInner=new AliExternalTrackParam(xmin,alpha,p,c); | |
468 | Quadratic2Helix(a,ca,b,cb,xmax,p,c); | |
469 | fOuter=new AliExternalTrackParam(xmax,alpha,p,c); | |
470 | } | |
471 | ||
472 | void AliTPCTracklet::Quadratic2Helix(Double_t *a,Double_t *ca, | |
473 | Double_t *b,Double_t *cb, | |
474 | Double_t x0, | |
475 | Double_t *p,Double_t *c) { | |
476 | // y(x)=a[0]+a[1]*x+a[2]*x^2 | |
477 | // z(x)=b[0]+b[1]*x | |
478 | // parametrises the corosponding helix at x0 | |
479 | ||
480 | // get the polynoms at x0 | |
481 | Double_t a0=x0*x0*a[2] + x0*a[1] + a[0]; | |
482 | Double_t a1=2.*x0*a[2] + a[1]; | |
483 | Double_t a2= a[2]; | |
484 | Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5]))); | |
485 | Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5])); | |
486 | Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]); | |
487 | Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]); | |
488 | Double_t ca21=ca[3]+x0*2.*ca[5]; | |
489 | Double_t ca22=ca[5]; | |
490 | ||
491 | Double_t b0=x0*b[1] + b[0]; | |
492 | Double_t b1= b[1]; | |
493 | Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]); | |
494 | Double_t cb10=cb[1]+x0*cb[2]; | |
495 | Double_t cb11=cb[2]; | |
496 | ||
497 | // transform to helix parameters | |
498 | Double_t f =1.+a1*a1; | |
499 | Double_t f2 =f*f; | |
500 | Double_t fi =1./f; | |
501 | Double_t fi12=TMath::Sqrt(fi); | |
502 | Double_t fi32=fi*fi12; | |
503 | Double_t fi2 =fi*fi; | |
504 | Double_t fi52=fi2*fi12; | |
505 | Double_t fi3 =fi2*fi; | |
506 | Double_t fi5 =fi2*fi3; | |
507 | ||
508 | Double_t xyz[3]={0.}; // TODO... | |
509 | Double_t fc=1./(GetBz(xyz)*kB2C); | |
510 | ||
511 | p[0]=a0; // y0 | |
512 | p[1]=b0; // z0 | |
513 | p[2]=a1*fi12; // snp | |
514 | p[3]=b1; // tgl | |
515 | p[4]=2.*a2*fi32*fc; // 1/pt | |
516 | ||
517 | c[0] =ca00; // y0-y0 | |
518 | c[1] =0.; // z0-y0 | |
519 | c[2] =cb00; // z0-z0 | |
520 | c[3] =ca10*fi32; // snp-y0 | |
521 | c[4] =0.; // snp-z0 | |
522 | c[5] =ca11*fi3; // snp-snp | |
523 | c[6] =0.; // tgl-y0 | |
524 | c[7] =cb10; // tgl-z0 | |
525 | c[8] =0.; // tgl-snp | |
526 | c[9] =cb11; // tgl-tgl | |
527 | c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc; // 1/pt-y0 | |
528 | c[11]=0.; // 1/pt-z0 | |
529 | c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp | |
530 | c[13]=0.; // 1/pt-tgl | |
531 | c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5 | |
532 | *fc*fc; // 1/pt-1/pt | |
533 | } | |
534 | ||
535 | ||
536 | void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) { | |
537 | TLinearFitter fy(2); | |
538 | fy.StoreData(kFALSE); | |
539 | fy.SetFormula("hyp2"); | |
540 | Double_t xmax=-1.; | |
541 | Double_t xmin=1000.; | |
542 | for (Int_t i=0;i<160;++i) { | |
543 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
967eae0d | 544 | if (c && RejectCluster(c)) continue; |
9318a5b4 | 545 | if (c&&c->GetDetector()==sector) { |
546 | Double_t x=c->GetX(); | |
547 | Double_t y=c->GetY(); | |
548 | Double_t xy[2]={x,y}; | |
549 | Double_t r=x*x+y*y; | |
550 | Double_t errx=1.,erry=1.;//TODO! | |
551 | Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry); | |
552 | err=1.; | |
553 | fy.AddPoint(xy,r,err); | |
554 | xmax=TMath::Max(xmax,x); | |
555 | xmin=TMath::Min(xmin,x); | |
556 | } | |
557 | } | |
558 | fy.Eval(); | |
559 | Double_t a[3]={fy.GetParameter(0), | |
560 | fy.GetParameter(1), | |
561 | fy.GetParameter(2)}; | |
562 | Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0), | |
563 | fy.GetCovarianceMatrixElement(1,0), | |
564 | fy.GetCovarianceMatrixElement(1,1), | |
565 | fy.GetCovarianceMatrixElement(2,0), | |
566 | fy.GetCovarianceMatrixElement(2,1), | |
567 | fy.GetCovarianceMatrixElement(2,2)}; | |
568 | ||
569 | TLinearFitter fz(1); | |
570 | fz.StoreData(kFALSE); | |
571 | fz.SetFormula("hyp1"); | |
572 | Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]); | |
573 | Double_t oldx=0.; | |
574 | Double_t oldy=R; | |
575 | Double_t phi=0.; | |
576 | for (Int_t i=0;i<160;++i) { | |
577 | AliTPCclusterMI *c=track->GetClusterPointer(i); | |
967eae0d | 578 | if (c && RejectCluster(c)) continue; |
9318a5b4 | 579 | if (c&&c->GetDetector()==sector) { |
580 | Double_t x=c->GetX(); | |
581 | Double_t y=c->GetY(); | |
582 | Double_t dx=x-oldx; | |
583 | Double_t dy=y-oldy; | |
584 | phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R)); | |
585 | Double_t err=1.; | |
586 | fz.AddPoint(&phi,c->GetZ(),err); | |
587 | oldx=x; | |
588 | oldy=y; | |
589 | } | |
590 | } | |
591 | fz.Eval(); | |
592 | Double_t b[2]={fz.GetParameter(0), | |
593 | fz.GetParameter(1)}; | |
594 | Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0), | |
595 | fz.GetCovarianceMatrixElement(1,0), | |
596 | fz.GetCovarianceMatrixElement(1,1)}; | |
597 | ||
598 | Double_t p[5]; | |
599 | Double_t c[15]; | |
600 | Double_t alpha=track->GetAlpha(); | |
601 | if (Riemann2Helix(a,ca,b,cb,0.,p,c)) | |
602 | fPrimary=new AliExternalTrackParam(0.,alpha,p,c); | |
603 | if (Riemann2Helix(a,ca,b,cb,xmin,p,c)) | |
604 | fInner=new AliExternalTrackParam(xmin,alpha,p,c); | |
605 | if (Riemann2Helix(a,ca,b,cb,xmax,p,c)) | |
606 | fOuter=new AliExternalTrackParam(xmax,alpha,p,c); | |
607 | } | |
608 | ||
c94a79e1 | 609 | Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/, |
610 | Double_t *b,Double_t */*cb*/, | |
9318a5b4 | 611 | Double_t x0, |
612 | Double_t *p,Double_t *c) { | |
613 | //TODO: signs! | |
614 | ||
615 | Double_t xr0=.5*a[1]; | |
616 | Double_t yr0=.5*a[2]; | |
617 | Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]); | |
618 | Double_t dx=x0-xr0; | |
619 | if (dx*dx>=R*R) return kFALSE; | |
60e55aee | 620 | Double_t dy=TMath::Sqrt((R-dx)*(R+dx)); //sign!! |
9318a5b4 | 621 | if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy)) |
622 | dy=-dy; | |
623 | Double_t y0=yr0+dy; | |
624 | Double_t tgp=-dx/dy; //TODO: dy!=0 | |
625 | Double_t z0=b[0]+TMath::ATan(tgp)*b[1]; | |
626 | Double_t xyz[3]={x0,y0,z0}; | |
627 | Double_t fc=1./(GetBz(xyz)*kB2C); | |
628 | fc=1; | |
629 | p[0]=y0; // y0 | |
630 | p[1]=z0; // z0 | |
631 | p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp | |
632 | p[3]=b[1]; // tgl | |
633 | p[4]=1./R*fc; // 1/pt | |
634 | ||
635 | c[0] =0.; // y0-y0 | |
636 | c[1] =0.; // z0-y0 | |
637 | c[2] =0.; // z0-z0 | |
638 | c[3] =0.; // snp-y0 | |
639 | c[4] =0.; // snp-z0 | |
640 | c[5] =0.; // snp-snp | |
641 | c[6] =0.; // tgl-y0 | |
642 | c[7] =0.; // tgl-z0 | |
643 | c[8] =0.; // tgl-snp | |
644 | c[9] =0.; // tgl-tgl | |
645 | c[10]=0.; // 1/pt-y0 | |
646 | c[11]=0.; // 1/pt-z0 | |
647 | c[12]=0.; // 1/pt-snp | |
648 | c[13]=0.; // 1/pt-tgl | |
649 | c[14]=0.; // 1/pt-1/pt | |
650 | ||
651 | return kTRUE; | |
652 | } | |
653 | ||
654 | TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track, | |
655 | TrackType type, | |
656 | Bool_t storeClusters, | |
69209e2c | 657 | Int_t minClusters, |
658 | Int_t maxTracklets) { | |
659 | // The tracklet factory: It creates several tracklets out of a track. They | |
660 | // are created for sectors that fullfill the constraint of having enough | |
661 | // clusters inside. Futhermore you can specify the maximum amount of | |
662 | // tracklets that are to be created. | |
663 | // The tracklets appear in a sorted fashion, beginning with those having the | |
664 | // most clusters. | |
665 | ||
666 | Int_t sectors[72]={0}; | |
667 | for (Int_t i=0;i<160;++i) { | |
9318a5b4 | 668 | AliTPCclusterMI *c=track->GetClusterPointer(i); |
967eae0d | 669 | if (c && RejectCluster(c)) continue; |
69209e2c | 670 | if (c) |
671 | ++sectors[c->GetDetector()]; | |
672 | } | |
673 | Int_t indices[72]; | |
674 | TMath::Sort(72,sectors,indices); | |
675 | TObjArray tracklets; | |
676 | if (maxTracklets>72) maxTracklets=72; // just to protect against "users". | |
677 | for (Int_t i=0;i<maxTracklets&§ors[indices[i]]>=minClusters;++i) { | |
9318a5b4 | 678 | tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters)); |
69209e2c | 679 | } |
680 | return tracklets; | |
681 | } | |
9318a5b4 | 682 | |
c94a79e1 | 683 | TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &/*clusters*/, |
684 | TrackType /*type*/, | |
685 | Bool_t /*storeClusters*/, | |
686 | Int_t /*minClusters*/, | |
687 | Int_t /*maxTracklets*/) { | |
688 | // TODO! | |
689 | ||
71d2416b | 690 | TObjArray tracklets; |
691 | return tracklets; | |
9318a5b4 | 692 | } |
693 | ||
694 | Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1, | |
695 | const AliTPCTracklet &t2, | |
696 | AliExternalTrackParam *&t1m, | |
697 | AliExternalTrackParam *&t2m) { | |
698 | // This function propagates two Tracklets to a common x-coordinate. This | |
699 | // x is dermined as the one that is in the middle of the two tracklets (they | |
700 | // are assumed to live on two distinct x-intervalls). | |
701 | // The inner parametrisation of the outer Tracklet and the outer | |
702 | // parametrisation of the inner Tracklet are used and propagated to this | |
703 | // common x. This result is saved not inside the Tracklets but two new | |
704 | // ExternalTrackParams are created (that means you might want to delete | |
705 | // them). | |
706 | // In the case that the alpha angles of the Tracklets differ both angles | |
707 | // are tried out for this propagation. | |
708 | // In case of any failure kFALSE is returned, no AliExternalTrackParam | |
709 | // is created und the pointers are set to 0. | |
710 | ||
711 | if (t1.GetInner() && t1.GetOuter() && | |
712 | t2.GetInner() && t2.GetOuter()) { | |
713 | if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) { | |
714 | t1m=new AliExternalTrackParam(*t1.GetOuter()); | |
715 | t2m=new AliExternalTrackParam(*t2.GetInner()); | |
716 | } | |
717 | else { | |
718 | t1m=new AliExternalTrackParam(*t1.GetInner()); | |
719 | t2m=new AliExternalTrackParam(*t2.GetOuter()); | |
720 | } | |
721 | Double_t mx=.5*(t1m->GetX()+t2m->GetX()); | |
273ccc8a | 722 | //Double_t b1,b2; |
9318a5b4 | 723 | Double_t xyz[3]; |
724 | t1m->GetXYZ(xyz); | |
273ccc8a | 725 | //b1=GetBz(xyz); |
726 | Double_t b1[3]; AliTracker::GetBxByBz(xyz,b1); | |
9318a5b4 | 727 | t2m->GetXYZ(xyz); |
273ccc8a | 728 | //b2=GetBz(xyz); |
729 | Double_t b2[3]; AliTracker::GetBxByBz(xyz,b2); | |
9318a5b4 | 730 | if (t1m->Rotate(t2m->GetAlpha()) |
273ccc8a | 731 | //&& t1m->PropagateTo(mx,b1) |
732 | //&& t2m->PropagateTo(mx,b2)); | |
733 | && t1m->PropagateToBxByBz(mx,b1) | |
734 | && t2m->PropagateToBxByBz(mx,b2)); | |
9318a5b4 | 735 | else |
736 | if (t2m->Rotate(t1m->GetAlpha()) | |
273ccc8a | 737 | //&& t1m->PropagateTo(mx,b1) |
738 | //&& t2m->PropagateTo(mx,b2)); | |
739 | && t1m->PropagateToBxByBz(mx,b1) | |
740 | && t2m->PropagateToBxByBz(mx,b2)); | |
9318a5b4 | 741 | else { |
742 | delete t1m; | |
743 | delete t2m; | |
744 | t1m=t2m=0; | |
745 | } | |
746 | } | |
747 | else { | |
748 | t1m=t2m=0; | |
749 | } | |
750 | return t1m&&t2m; | |
751 | } | |
752 | ||
f7a1cc68 | 753 | double AliTPCTracklet::GetBz(Double_t *xyz) |
754 | { | |
755 | return AliTracker::GetBz(xyz); | |
9318a5b4 | 756 | } |
757 | ||
758 | void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c, | |
759 | Double_t *x) { | |
760 | // This function generates a n-dimensional random variable x with mean | |
761 | // p and covariance c. | |
762 | // That is done using the cholesky decomposition of the covariance matrix, | |
763 | // Begin_Latex C=U^{t} U End_Latex, with Begin_Latex U End_Latex being an | |
764 | // upper triangular matrix. Given a vector v of iid gausian random variables | |
765 | // with variance 1 one obtains the asked result as: Begin_Latex x=U^t v | |
766 | // End_Latex. | |
767 | // c is expected to be in a lower triangular format: | |
768 | // c[0] | |
769 | // c[1] c[2] | |
770 | // c[3] c[4] c[5] | |
771 | // etc. | |
772 | static TRandom3 random; | |
773 | Double_t *c2= new Double_t[ndim*ndim]; | |
774 | Int_t k=0; | |
775 | for (Int_t i=0;i<ndim;++i) | |
776 | for (Int_t j=0;j<=i;++j) | |
777 | c2[i*ndim+j]=c2[j*ndim+i]=c[k++]; | |
778 | TMatrixDSym cm(ndim,c2); | |
779 | delete[] c2; | |
780 | TDecompChol chol(cm); | |
781 | chol.Decompose(); | |
782 | const TVectorD pv(ndim); | |
783 | const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p)); | |
784 | TVectorD xv(ndim); | |
785 | xv.Use(ndim,x); | |
786 | for (Int_t i=0;i<ndim;++i) | |
787 | xv[i]=random.Gaus(); | |
788 | TMatrixD L=chol.GetU(); | |
789 | L.T(); | |
790 | xv=L*xv+pv; | |
791 | } | |
792 | ||
793 | TEllipse AliTPCTracklet::ErrorEllipse(Double_t x,Double_t y, | |
794 | Double_t sx,Double_t sy,Double_t sxy) { | |
795 | /* Begin_Latex | |
796 | r_{1,2}=1/2 (a+c#pm#sqrt{(a-c)^{2}+(2b)^{2}}) | |
797 | End_Latex */ | |
798 | Double_t det1=1./(sx*sy-sxy*sxy); | |
799 | Double_t a=sy*det1; | |
800 | Double_t b=-sxy*det1; | |
801 | Double_t c=sx*det1; | |
802 | Double_t d=c-a; | |
803 | Double_t s=TMath::Sqrt(d*d+4.*b*b); | |
804 | Double_t r1=TMath::Sqrt(.5*(a+c-s)); | |
805 | Double_t r2=TMath::Sqrt(.5*(a+c+s)); | |
806 | Double_t alpha=.5*TMath::ATan2(2.*b,d); | |
807 | return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg()); | |
808 | } | |
809 | ||
810 | void AliTPCTracklet::Test(const char* filename) { | |
811 | /* | |
812 | aliroot | |
813 | AliTPCTracklet::Test(""); | |
814 | TFile f("AliTPCTrackletDebug.root"); | |
815 | TTree *t=f.Get("AliTPCTrackletDebug"); | |
816 | t->Draw("p0:p4"); | |
817 | TEllipse e=AliTPCTracklet::ErrorEllipse(0.,0.,4.,1.,1.8); | |
818 | e.Draw(); | |
819 | */ | |
c94a79e1 | 820 | TTreeSRedirector ds(filename); |
9318a5b4 | 821 | Double_t p[5]={0.}; |
822 | Double_t c[15]={4., | |
823 | 0.,4., | |
824 | 0.,0.,9., | |
825 | 0.,0.,0.,16., | |
826 | 1.8,0.,0.,0.,1.}; | |
827 | for (Int_t i=0;i<10000;++i) { | |
828 | Double_t x[5]; | |
829 | RandomND(5,p,c,x); | |
830 | ds<<"AliTPCTrackletDebug" | |
831 | <<"p0="<<x[0] | |
832 | <<"p1="<<x[1] | |
833 | <<"p2="<<x[2] | |
834 | <<"p3="<<x[3] | |
835 | <<"p4="<<x[4] | |
836 | <<"\n"; | |
837 | } | |
838 | ||
839 | /* | |
840 | Double_t b; | |
841 | Double_t x=0.; | |
842 | Double_t alpha=0.; | |
843 | Double_t param[5]={0.}; | |
844 | Double_t covar[15]={1., | |
845 | 0.,1., | |
846 | 0.,0.,1., | |
847 | 0.,0.,0.,1., | |
848 | 0.,0.,0.,0.,1.}; | |
849 | AliExternalTrackParam track(x,alpha,param,covar); | |
850 | ||
851 | ||
852 | ||
853 | for (Int_t i=0;i<points.GetNPoints();++i) { | |
854 | Double_t x=0.; | |
855 | Double_t alpha=0.; | |
856 | Double_t param[5]={0.}; | |
857 | Double_t covar[15]={1., | |
858 | 0.,1., | |
859 | 0.,0.,1., | |
860 | 0.,0.,0.,1., | |
861 | 0.,0.,0.,0.,1.}; | |
862 | AliExternalTrackParam track(x,alpha,param,covar); | |
863 | for (x=90.;x<250.;x+=1.) { | |
864 | track.PropagateTo(x,b); | |
865 | AliTPCclusterMI c(); | |
866 | } | |
867 | } | |
868 | */ | |
869 | } | |
967eae0d | 870 | |
871 | ||
872 | Bool_t AliTPCTracklet::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){ | |
873 | // | |
874 | // check the acceptance of cluster | |
875 | // Cut on edge effects | |
876 | // | |
877 | Bool_t isReject = kFALSE; | |
878 | Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18); | |
879 | Float_t dist = edgeY - TMath::Abs(cl->GetY()); | |
880 | if (param) dist = edgeY - TMath::Abs(param->GetY()); | |
881 | if (dist<fgEdgeCutY) isReject=kTRUE; | |
882 | if (cl->GetType()<0) isReject=kTRUE; | |
883 | return isReject; | |
884 | } |