]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCTracklet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
coverity fix
[u/mrichter/AliRoot.git] / TPC / AliTPCTracklet.cxx
CommitLineData
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>
35using namespace std;
69209e2c 36
37ClassImp(AliTPCTracklet)
38
9318a5b4 39const Double_t AliTPCTracklet::kB2C=0.299792458e-3;
967eae0d 40Float_t AliTPCTracklet::fgEdgeCutY=3;
41Float_t AliTPCTracklet::fgEdgeCutX=0;
9318a5b4 42
69209e2c 43AliTPCTracklet::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 51AliTPCTracklet::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 95AliTPCTracklet::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
102AliTPCTracklet::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
123AliTPCTracklet& 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;
9318a5b4 129 delete fClusters;
130 if (t.fClusters) {
131 fClusters=new AliTPCclusterMI[t.fNStoredClusters];
132 for (int i=0;i<t.fNStoredClusters;++i)
133 fClusters[i]=t.fClusters[i];
134 }
135 else
136 fClusters=0;
137 fSector=t.fSector;
69209e2c 138 if (t.fOuter) {
139 if (fOuter)
140 *fOuter=*t.fOuter;
141 else
142 fOuter=new AliExternalTrackParam(*t.fOuter);
143 }
144 else {
145 delete fOuter;
146 fOuter=0;
147 }
148
149 if (t.fInner) {
150 if (fInner)
151 *fInner=*t.fInner;
152 else
153 fInner=new AliExternalTrackParam(*t.fInner);
154 }
155 else {
156 delete fInner;
157 fInner=0;
158 }
159
160 if (t.fPrimary) {
161 if (fPrimary)
162 *fPrimary=*t.fPrimary;
163 else
164 fPrimary=new AliExternalTrackParam(*t.fPrimary);
165 }
166 else {
167 delete fPrimary;
168 fPrimary=0;
169 }
170 }
171 return *this;
172}
173
174AliTPCTracklet::~AliTPCTracklet() {
175 //
176 // The destructor. Yes, you can even destruct tracklets.
177 //
9318a5b4 178 delete fClusters;
69209e2c 179 delete fOuter;
180 delete fInner;
181 delete fPrimary;
182}
183
ef0a7e89 184
185
ef0a7e89 186
187
188void AliTPCTracklet::FitKalman(const AliTPCseed *seed,Int_t sector) {
189 //
190 // Fit using Kalman filter
191 //
192 AliTPCseed *track=new AliTPCseed(*seed);
193 if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
194 delete track;
195 return;
196 }
197 // fit from inner to outer row
198 Double_t covar[15];
199 for (Int_t i=0;i<15;i++) covar[i]=0;
b322e06a 200 covar[0]=1.*1.;
201 covar[2]=1.*1.;
202 covar[5]=1.*1./(64.*64.);
203 covar[9]=1.*1./(64.*64.);
204 covar[14]=0; // keep pt
ef0a7e89 205 Float_t xmin=1000, xmax=-10000;
206 Int_t imin=158, imax=0;
207 for (Int_t i=0;i<160;i++) {
208 AliTPCclusterMI *c=track->GetClusterPointer(i);
209 if (!c) continue;
210 if (c->GetDetector()!=sector) continue;
211 if (c->GetX()<xmin) xmin=c->GetX();
212 if (c->GetX()>xmax) xmax=c->GetX();
213 if (i<imin) imin=i;
214 if (i>imax) imax=i;
215 }
216 if(imax-imin<10) {
217 delete track;
218 return;
219 }
220
221 for (Float_t x=track->GetX(); x<xmin; x++) track->PropagateTo(x);
222 track->AddCovariance(covar);
223 //
224 AliExternalTrackParam paramIn;
225 AliExternalTrackParam paramOut;
226 Bool_t isOK=kTRUE;
227 //
228 //
229 //
230 for (Int_t i=imin; i<=imax; i++){
231 AliTPCclusterMI *c=track->GetClusterPointer(i);
232 if (!c) continue;
ef0a7e89 233 Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
234 Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
91350f5b 235 AliTPCseed::GetError(c, track,cov[0],cov[2]);
236 cov[0]*=cov[0];
237 cov[2]*=cov[2];
ef0a7e89 238 if (!track->PropagateTo(r[0])) {
239 isOK=kFALSE;
240 break;
241 }
42b40d07 242 if (RejectCluster(c,track)) continue;
ef0a7e89 243 if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
244 }
245 if (!isOK) { delete track; return;}
246 track->AddCovariance(covar);
247 //
248 //
249 //
250 for (Int_t i=imax; i>=imin; i--){
251 AliTPCclusterMI *c=track->GetClusterPointer(i);
252 if (!c) continue;
ef0a7e89 253 Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
91350f5b 254 Double_t cov[3]={0.01,0.,0.01};
255 AliTPCseed::GetError(c, track,cov[0],cov[2]);
256 cov[0]*=cov[0];
257 cov[2]*=cov[2];
ef0a7e89 258 if (!track->PropagateTo(r[0])) {
259 isOK=kFALSE;
260 break;
261 }
42b40d07 262 if (RejectCluster(c,track)) continue;
ef0a7e89 263 if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
264 }
265 if (!isOK) { delete track; return;}
266 paramIn = *track;
267 track->AddCovariance(covar);
268 //
269 //
270 for (Int_t i=imin; i<=imax; i++){
271 AliTPCclusterMI *c=track->GetClusterPointer(i);
272 if (!c) continue;
ef0a7e89 273 Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
91350f5b 274 Double_t cov[3]={0.01,0.,0.01};
275 AliTPCseed::GetError(c, track,cov[0],cov[2]);
276 cov[0]*=cov[0];
277 cov[2]*=cov[2];
ef0a7e89 278 if (!track->PropagateTo(r[0])) {
279 isOK=kFALSE;
280 break;
281 }
42b40d07 282 if (RejectCluster(c,track)) continue;
ef0a7e89 283 if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
284 }
285 if (!isOK) { delete track; return;}
286 paramOut=*track;
287 //
288 //
289 //
290 fOuter=new AliExternalTrackParam(paramOut);
291 fInner=new AliExternalTrackParam(paramIn);
292 //
ef0a7e89 293 delete track;
294}
295
296
297
9318a5b4 298
299void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
300 TrackType type) {
301 TLinearFitter fy(1);
302 TLinearFitter fz(1);
303 fy.StoreData(kFALSE);
304 fz.StoreData(kFALSE);
305 switch (type) {
306 case kLinear:
307 fy.SetFormula("1 ++ x");
308 fz.SetFormula("1 ++ x");
309 break;
310 case kQuadratic:
311 fy.SetFormula("1 ++ x ++ x*x");
312 fz.SetFormula("1 ++ x");
313 break;
c94a79e1 314 case kKalman:
315 case kRiemann:
316 break;
9318a5b4 317 }
318 Double_t xmax=-1.;
319 Double_t xmin=1000.;
320 for (Int_t i=0;i<160;++i) {
321 AliTPCclusterMI *c=track->GetClusterPointer(i);
967eae0d 322 if (c && RejectCluster(c)) continue;
9318a5b4 323 if (c&&c->GetDetector()==sector) {
324 Double_t x=c->GetX();
325 fy.AddPoint(&x,c->GetY());
326 fz.AddPoint(&x,c->GetZ());
327 xmax=TMath::Max(xmax,x);
328 xmin=TMath::Min(xmin,x);
329 }
330 }
331 fy.Eval();
332 fz.Eval();
333 Double_t a[3]={fy.GetParameter(0),
334 fy.GetParameter(1),
335 type==kQuadratic?fy.GetParameter(2):0.};
336 Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
337 fy.GetCovarianceMatrixElement(1,0),
338 fy.GetCovarianceMatrixElement(1,1),
339 type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0.,
340 type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0.,
341 type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.};
342 for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters;
343 Double_t b[2]={fz.GetParameter(0),
344 fz.GetParameter(1)};
345 Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
346 fz.GetCovarianceMatrixElement(1,0),
347 fz.GetCovarianceMatrixElement(1,1)};
348 for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters;
349 Double_t p[5];
350 Double_t c[15];
351 Double_t alpha=track->GetAlpha();
352 Quadratic2Helix(a,ca,b,cb,0.,p,c);
353 fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
354 Quadratic2Helix(a,ca,b,cb,xmin,p,c);
355 fInner=new AliExternalTrackParam(xmin,alpha,p,c);
356 Quadratic2Helix(a,ca,b,cb,xmax,p,c);
357 fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
358}
359
360void AliTPCTracklet::Quadratic2Helix(Double_t *a,Double_t *ca,
361 Double_t *b,Double_t *cb,
362 Double_t x0,
363 Double_t *p,Double_t *c) {
364 // y(x)=a[0]+a[1]*x+a[2]*x^2
365 // z(x)=b[0]+b[1]*x
366 // parametrises the corosponding helix at x0
367
368 // get the polynoms at x0
369 Double_t a0=x0*x0*a[2] + x0*a[1] + a[0];
370 Double_t a1=2.*x0*a[2] + a[1];
371 Double_t a2= a[2];
372 Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5])));
373 Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5]));
374 Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]);
375 Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]);
376 Double_t ca21=ca[3]+x0*2.*ca[5];
377 Double_t ca22=ca[5];
378
379 Double_t b0=x0*b[1] + b[0];
380 Double_t b1= b[1];
381 Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]);
382 Double_t cb10=cb[1]+x0*cb[2];
383 Double_t cb11=cb[2];
384
385 // transform to helix parameters
386 Double_t f =1.+a1*a1;
387 Double_t f2 =f*f;
388 Double_t fi =1./f;
389 Double_t fi12=TMath::Sqrt(fi);
390 Double_t fi32=fi*fi12;
391 Double_t fi2 =fi*fi;
392 Double_t fi52=fi2*fi12;
393 Double_t fi3 =fi2*fi;
394 Double_t fi5 =fi2*fi3;
395
396 Double_t xyz[3]={0.}; // TODO...
397 Double_t fc=1./(GetBz(xyz)*kB2C);
398
399 p[0]=a0; // y0
400 p[1]=b0; // z0
401 p[2]=a1*fi12; // snp
402 p[3]=b1; // tgl
403 p[4]=2.*a2*fi32*fc; // 1/pt
404
405 c[0] =ca00; // y0-y0
406 c[1] =0.; // z0-y0
407 c[2] =cb00; // z0-z0
408 c[3] =ca10*fi32; // snp-y0
409 c[4] =0.; // snp-z0
410 c[5] =ca11*fi3; // snp-snp
411 c[6] =0.; // tgl-y0
412 c[7] =cb10; // tgl-z0
413 c[8] =0.; // tgl-snp
414 c[9] =cb11; // tgl-tgl
415 c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc; // 1/pt-y0
416 c[11]=0.; // 1/pt-z0
417 c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp
418 c[13]=0.; // 1/pt-tgl
419 c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5
420 *fc*fc; // 1/pt-1/pt
421}
422
423
424void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) {
425 TLinearFitter fy(2);
426 fy.StoreData(kFALSE);
427 fy.SetFormula("hyp2");
428 Double_t xmax=-1.;
429 Double_t xmin=1000.;
430 for (Int_t i=0;i<160;++i) {
431 AliTPCclusterMI *c=track->GetClusterPointer(i);
967eae0d 432 if (c && RejectCluster(c)) continue;
9318a5b4 433 if (c&&c->GetDetector()==sector) {
434 Double_t x=c->GetX();
435 Double_t y=c->GetY();
436 Double_t xy[2]={x,y};
437 Double_t r=x*x+y*y;
438 Double_t errx=1.,erry=1.;//TODO!
439 Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry);
440 err=1.;
441 fy.AddPoint(xy,r,err);
442 xmax=TMath::Max(xmax,x);
443 xmin=TMath::Min(xmin,x);
444 }
445 }
446 fy.Eval();
447 Double_t a[3]={fy.GetParameter(0),
448 fy.GetParameter(1),
449 fy.GetParameter(2)};
450 Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
451 fy.GetCovarianceMatrixElement(1,0),
452 fy.GetCovarianceMatrixElement(1,1),
453 fy.GetCovarianceMatrixElement(2,0),
454 fy.GetCovarianceMatrixElement(2,1),
455 fy.GetCovarianceMatrixElement(2,2)};
456
457 TLinearFitter fz(1);
458 fz.StoreData(kFALSE);
459 fz.SetFormula("hyp1");
460 Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
461 Double_t oldx=0.;
462 Double_t oldy=R;
463 Double_t phi=0.;
464 for (Int_t i=0;i<160;++i) {
465 AliTPCclusterMI *c=track->GetClusterPointer(i);
967eae0d 466 if (c && RejectCluster(c)) continue;
9318a5b4 467 if (c&&c->GetDetector()==sector) {
468 Double_t x=c->GetX();
469 Double_t y=c->GetY();
470 Double_t dx=x-oldx;
471 Double_t dy=y-oldy;
472 phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R));
473 Double_t err=1.;
474 fz.AddPoint(&phi,c->GetZ(),err);
475 oldx=x;
476 oldy=y;
477 }
478 }
479 fz.Eval();
480 Double_t b[2]={fz.GetParameter(0),
481 fz.GetParameter(1)};
482 Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
483 fz.GetCovarianceMatrixElement(1,0),
484 fz.GetCovarianceMatrixElement(1,1)};
485
486 Double_t p[5];
487 Double_t c[15];
488 Double_t alpha=track->GetAlpha();
489 if (Riemann2Helix(a,ca,b,cb,0.,p,c))
490 fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
491 if (Riemann2Helix(a,ca,b,cb,xmin,p,c))
492 fInner=new AliExternalTrackParam(xmin,alpha,p,c);
493 if (Riemann2Helix(a,ca,b,cb,xmax,p,c))
494 fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
495}
496
c94a79e1 497Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/,
498 Double_t *b,Double_t */*cb*/,
9318a5b4 499 Double_t x0,
500 Double_t *p,Double_t *c) {
501 //TODO: signs!
502
503 Double_t xr0=.5*a[1];
504 Double_t yr0=.5*a[2];
505 Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
506 Double_t dx=x0-xr0;
507 if (dx*dx>=R*R) return kFALSE;
60e55aee 508 Double_t dy=TMath::Sqrt((R-dx)*(R+dx)); //sign!!
9318a5b4 509 if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
510 dy=-dy;
511 Double_t y0=yr0+dy;
512 Double_t tgp=-dx/dy; //TODO: dy!=0
513 Double_t z0=b[0]+TMath::ATan(tgp)*b[1];
514 Double_t xyz[3]={x0,y0,z0};
515 Double_t fc=1./(GetBz(xyz)*kB2C);
516 fc=1;
517 p[0]=y0; // y0
518 p[1]=z0; // z0
519 p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp
520 p[3]=b[1]; // tgl
521 p[4]=1./R*fc; // 1/pt
522
523 c[0] =0.; // y0-y0
524 c[1] =0.; // z0-y0
525 c[2] =0.; // z0-z0
526 c[3] =0.; // snp-y0
527 c[4] =0.; // snp-z0
528 c[5] =0.; // snp-snp
529 c[6] =0.; // tgl-y0
530 c[7] =0.; // tgl-z0
531 c[8] =0.; // tgl-snp
532 c[9] =0.; // tgl-tgl
533 c[10]=0.; // 1/pt-y0
534 c[11]=0.; // 1/pt-z0
535 c[12]=0.; // 1/pt-snp
536 c[13]=0.; // 1/pt-tgl
537 c[14]=0.; // 1/pt-1/pt
538
539 return kTRUE;
540}
541
542TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
543 TrackType type,
544 Bool_t storeClusters,
69209e2c 545 Int_t minClusters,
546 Int_t maxTracklets) {
547// The tracklet factory: It creates several tracklets out of a track. They
548// are created for sectors that fullfill the constraint of having enough
549// clusters inside. Futhermore you can specify the maximum amount of
550// tracklets that are to be created.
551// The tracklets appear in a sorted fashion, beginning with those having the
552// most clusters.
553
554 Int_t sectors[72]={0};
555 for (Int_t i=0;i<160;++i) {
9318a5b4 556 AliTPCclusterMI *c=track->GetClusterPointer(i);
967eae0d 557 if (c && RejectCluster(c)) continue;
69209e2c 558 if (c)
559 ++sectors[c->GetDetector()];
560 }
561 Int_t indices[72];
562 TMath::Sort(72,sectors,indices);
563 TObjArray tracklets;
564 if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
565 for (Int_t i=0;i<maxTracklets&&sectors[indices[i]]>=minClusters;++i) {
9318a5b4 566 tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
69209e2c 567 }
568 return tracklets;
569}
9318a5b4 570
c94a79e1 571TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &/*clusters*/,
572 TrackType /*type*/,
573 Bool_t /*storeClusters*/,
574 Int_t /*minClusters*/,
575 Int_t /*maxTracklets*/) {
576 // TODO!
577
71d2416b 578 TObjArray tracklets;
579 return tracklets;
9318a5b4 580}
581
582Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
583 const AliTPCTracklet &t2,
584 AliExternalTrackParam *&t1m,
585 AliExternalTrackParam *&t2m) {
586 // This function propagates two Tracklets to a common x-coordinate. This
587 // x is dermined as the one that is in the middle of the two tracklets (they
588 // are assumed to live on two distinct x-intervalls).
589 // The inner parametrisation of the outer Tracklet and the outer
590 // parametrisation of the inner Tracklet are used and propagated to this
591 // common x. This result is saved not inside the Tracklets but two new
592 // ExternalTrackParams are created (that means you might want to delete
593 // them).
594 // In the case that the alpha angles of the Tracklets differ both angles
595 // are tried out for this propagation.
596 // In case of any failure kFALSE is returned, no AliExternalTrackParam
597 // is created und the pointers are set to 0.
598
599 if (t1.GetInner() && t1.GetOuter() &&
600 t2.GetInner() && t2.GetOuter()) {
601 if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) {
602 t1m=new AliExternalTrackParam(*t1.GetOuter());
603 t2m=new AliExternalTrackParam(*t2.GetInner());
604 }
605 else {
606 t1m=new AliExternalTrackParam(*t1.GetInner());
607 t2m=new AliExternalTrackParam(*t2.GetOuter());
608 }
609 Double_t mx=.5*(t1m->GetX()+t2m->GetX());
273ccc8a 610 //Double_t b1,b2;
9318a5b4 611 Double_t xyz[3];
612 t1m->GetXYZ(xyz);
273ccc8a 613 //b1=GetBz(xyz);
614 Double_t b1[3]; AliTracker::GetBxByBz(xyz,b1);
9318a5b4 615 t2m->GetXYZ(xyz);
273ccc8a 616 //b2=GetBz(xyz);
617 Double_t b2[3]; AliTracker::GetBxByBz(xyz,b2);
9318a5b4 618 if (t1m->Rotate(t2m->GetAlpha())
273ccc8a 619 //&& t1m->PropagateTo(mx,b1)
620 //&& t2m->PropagateTo(mx,b2));
621 && t1m->PropagateToBxByBz(mx,b1)
622 && t2m->PropagateToBxByBz(mx,b2));
9318a5b4 623 else
624 if (t2m->Rotate(t1m->GetAlpha())
273ccc8a 625 //&& t1m->PropagateTo(mx,b1)
626 //&& t2m->PropagateTo(mx,b2));
627 && t1m->PropagateToBxByBz(mx,b1)
628 && t2m->PropagateToBxByBz(mx,b2));
9318a5b4 629 else {
630 delete t1m;
631 delete t2m;
632 t1m=t2m=0;
633 }
634 }
635 else {
636 t1m=t2m=0;
637 }
638 return t1m&&t2m;
639}
640
f7a1cc68 641double AliTPCTracklet::GetBz(Double_t *xyz)
642{
643 return AliTracker::GetBz(xyz);
9318a5b4 644}
645
646void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c,
647 Double_t *x) {
648 // This function generates a n-dimensional random variable x with mean
649 // p and covariance c.
650 // That is done using the cholesky decomposition of the covariance matrix,
651 // Begin_Latex C=U^{t} U End_Latex, with Begin_Latex U End_Latex being an
652 // upper triangular matrix. Given a vector v of iid gausian random variables
653 // with variance 1 one obtains the asked result as: Begin_Latex x=U^t v
654 // End_Latex.
655 // c is expected to be in a lower triangular format:
656 // c[0]
657 // c[1] c[2]
658 // c[3] c[4] c[5]
659 // etc.
660 static TRandom3 random;
661 Double_t *c2= new Double_t[ndim*ndim];
662 Int_t k=0;
663 for (Int_t i=0;i<ndim;++i)
664 for (Int_t j=0;j<=i;++j)
665 c2[i*ndim+j]=c2[j*ndim+i]=c[k++];
666 TMatrixDSym cm(ndim,c2);
667 delete[] c2;
668 TDecompChol chol(cm);
669 chol.Decompose();
670 const TVectorD pv(ndim);
671 const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p));
672 TVectorD xv(ndim);
673 xv.Use(ndim,x);
674 for (Int_t i=0;i<ndim;++i)
675 xv[i]=random.Gaus();
676 TMatrixD L=chol.GetU();
677 L.T();
678 xv=L*xv+pv;
679}
680
681TEllipse AliTPCTracklet::ErrorEllipse(Double_t x,Double_t y,
682 Double_t sx,Double_t sy,Double_t sxy) {
683 /* Begin_Latex
684 r_{1,2}=1/2 (a+c#pm#sqrt{(a-c)^{2}+(2b)^{2}})
685 End_Latex */
686 Double_t det1=1./(sx*sy-sxy*sxy);
687 Double_t a=sy*det1;
688 Double_t b=-sxy*det1;
689 Double_t c=sx*det1;
690 Double_t d=c-a;
691 Double_t s=TMath::Sqrt(d*d+4.*b*b);
692 Double_t r1=TMath::Sqrt(.5*(a+c-s));
693 Double_t r2=TMath::Sqrt(.5*(a+c+s));
694 Double_t alpha=.5*TMath::ATan2(2.*b,d);
695 return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg());
696}
697
698void AliTPCTracklet::Test(const char* filename) {
699 /*
700 aliroot
701 AliTPCTracklet::Test("");
702 TFile f("AliTPCTrackletDebug.root");
703 TTree *t=f.Get("AliTPCTrackletDebug");
704 t->Draw("p0:p4");
705 TEllipse e=AliTPCTracklet::ErrorEllipse(0.,0.,4.,1.,1.8);
706 e.Draw();
707 */
c94a79e1 708 TTreeSRedirector ds(filename);
9318a5b4 709 Double_t p[5]={0.};
710 Double_t c[15]={4.,
711 0.,4.,
712 0.,0.,9.,
713 0.,0.,0.,16.,
714 1.8,0.,0.,0.,1.};
715 for (Int_t i=0;i<10000;++i) {
716 Double_t x[5];
717 RandomND(5,p,c,x);
718 ds<<"AliTPCTrackletDebug"
719 <<"p0="<<x[0]
720 <<"p1="<<x[1]
721 <<"p2="<<x[2]
722 <<"p3="<<x[3]
723 <<"p4="<<x[4]
724 <<"\n";
725 }
726
727 /*
728 Double_t b;
729 Double_t x=0.;
730 Double_t alpha=0.;
731 Double_t param[5]={0.};
732 Double_t covar[15]={1.,
733 0.,1.,
734 0.,0.,1.,
735 0.,0.,0.,1.,
736 0.,0.,0.,0.,1.};
737 AliExternalTrackParam track(x,alpha,param,covar);
738
739
740
741 for (Int_t i=0;i<points.GetNPoints();++i) {
742 Double_t x=0.;
743 Double_t alpha=0.;
744 Double_t param[5]={0.};
745 Double_t covar[15]={1.,
746 0.,1.,
747 0.,0.,1.,
748 0.,0.,0.,1.,
749 0.,0.,0.,0.,1.};
750 AliExternalTrackParam track(x,alpha,param,covar);
751 for (x=90.;x<250.;x+=1.) {
752 track.PropagateTo(x,b);
753 AliTPCclusterMI c();
754 }
755 }
756 */
757}
967eae0d 758
759
760Bool_t AliTPCTracklet::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
761 //
762 // check the acceptance of cluster
763 // Cut on edge effects
764 //
765 Bool_t isReject = kFALSE;
766 Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
767 Float_t dist = edgeY - TMath::Abs(cl->GetY());
768 if (param) dist = edgeY - TMath::Abs(param->GetY());
769 if (dist<fgEdgeCutY) isReject=kTRUE;
770 if (cl->GetType()<0) isReject=kTRUE;
771 return isReject;
772}