1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
24 #include "AliTPCTracklet.h"
25 #include "TObjArray.h"
26 #include "TLinearFitter.h"
27 #include "AliTPCseed.h"
28 #include "AliESDVertex.h"
29 #include "AliTracker.h"
30 #include "TTreeStream.h"
32 #include "TDecompChol.h"
37 ClassImp(AliTPCTracklet)
39 const Double_t AliTPCTracklet::kB2C=0.299792458e-3;
40 Float_t AliTPCTracklet::fgEdgeCutY=3;
41 Float_t AliTPCTracklet::fgEdgeCutX=0;
43 AliTPCTracklet::AliTPCTracklet()
44 : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0),
45 fInner(0),fPrimary(0) {
47 // The default constructor. It is intended to be used for I/O only.
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) {
56 // Contructor for a tracklet out of a track. Only clusters within a given
60 //TODO: only kalman works
62 for (Int_t i=0;i<160;++i) {
63 AliTPCclusterMI *c=track->GetClusterPointer(i);
64 if (c && RejectCluster(c)) continue;
65 if (c&&c->GetDetector()==sector)
70 fClusters=new AliTPCclusterMI[fNClusters];
71 for (Int_t i=0;i<160;++i) {
72 AliTPCclusterMI *c=track->GetClusterPointer(i);
73 if (c && RejectCluster(c)) continue;
74 if (c&&c->GetDetector()==sector)
75 fClusters[fNStoredClusters]=*c;
82 FitKalman(track,sector);
86 FitLinear(track,sector,type);
89 FitRiemann(track,sector);
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) {
102 AliTPCTracklet::AliTPCTracklet(const AliTPCTracklet &t)
103 : TObject(t),fNClusters(t.fNClusters),fNStoredClusters(t.fNStoredClusters),fClusters(0),
104 fSector(t.fSector),fOuter(0),fInner(0),
107 // The copy constructor. You can copy tracklets!
111 fClusters=new AliTPCclusterMI[t.fNStoredClusters];
112 for (int i=0;i<t.fNStoredClusters;++i)
113 fClusters[i]=t.fClusters[i];
116 fOuter=new AliExternalTrackParam(*t.fOuter);
118 fInner=new AliExternalTrackParam(*t.fInner);
120 fPrimary=new AliExternalTrackParam(*t.fPrimary);
123 AliTPCTracklet& AliTPCTracklet::operator=(const AliTPCTracklet &t) {
125 // The assignment constructor. You can assign tracklets!
128 fNClusters=t.fNClusters;
129 fNStoredClusters=fNStoredClusters;
132 fClusters=new AliTPCclusterMI[t.fNStoredClusters];
133 for (int i=0;i<t.fNStoredClusters;++i)
134 fClusters[i]=t.fClusters[i];
143 fOuter=new AliExternalTrackParam(*t.fOuter);
154 fInner=new AliExternalTrackParam(*t.fInner);
163 *fPrimary=*t.fPrimary;
165 fPrimary=new AliExternalTrackParam(*t.fPrimary);
175 AliTPCTracklet::~AliTPCTracklet() {
177 // The destructor. Yes, you can even destruct tracklets.
189 void AliTPCTracklet::FitKalman(const AliTPCseed *seed,Int_t sector) {
191 // Fit using Kalman filter
193 AliTPCseed *track=new AliTPCseed(*seed);
194 if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
198 // fit from inner to outer row
200 for (Int_t i=0;i<15;i++) covar[i]=0;
203 covar[5]=1.*1./(64.*64.);
204 covar[9]=1.*1./(64.*64.);
205 covar[14]=0; // keep pt
206 Float_t xmin=1000, xmax=-10000;
207 Int_t imin=158, imax=0;
208 for (Int_t i=0;i<160;i++) {
209 AliTPCclusterMI *c=track->GetClusterPointer(i);
211 if (c->GetDetector()!=sector) continue;
212 if (c->GetX()<xmin) xmin=c->GetX();
213 if (c->GetX()>xmax) xmax=c->GetX();
222 for (Float_t x=track->GetX(); x<xmin; x++) track->PropagateTo(x);
223 track->AddCovariance(covar);
225 AliExternalTrackParam paramIn;
226 AliExternalTrackParam paramOut;
231 for (Int_t i=imin; i<=imax; i++){
232 AliTPCclusterMI *c=track->GetClusterPointer(i);
234 Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
235 Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
236 AliTPCseed::GetError(c, track,cov[0],cov[2]);
239 if (!track->PropagateTo(r[0])) {
243 if (RejectCluster(c,track)) continue;
244 if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
246 if (!isOK) { delete track; return;}
247 track->AddCovariance(covar);
251 for (Int_t i=imax; i>=imin; i--){
252 AliTPCclusterMI *c=track->GetClusterPointer(i);
254 Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
255 Double_t cov[3]={0.01,0.,0.01};
256 AliTPCseed::GetError(c, track,cov[0],cov[2]);
259 if (!track->PropagateTo(r[0])) {
263 if (RejectCluster(c,track)) continue;
264 if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
266 if (!isOK) { delete track; return;}
268 track->AddCovariance(covar);
271 for (Int_t i=imin; i<=imax; i++){
272 AliTPCclusterMI *c=track->GetClusterPointer(i);
274 Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
275 Double_t cov[3]={0.01,0.,0.01};
276 AliTPCseed::GetError(c, track,cov[0],cov[2]);
279 if (!track->PropagateTo(r[0])) {
283 if (RejectCluster(c,track)) continue;
284 if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
286 if (!isOK) { delete track; return;}
291 fOuter=new AliExternalTrackParam(paramOut);
292 fInner=new AliExternalTrackParam(paramIn);
300 void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
304 fy.StoreData(kFALSE);
305 fz.StoreData(kFALSE);
308 fy.SetFormula("1 ++ x");
309 fz.SetFormula("1 ++ x");
312 fy.SetFormula("1 ++ x ++ x*x");
313 fz.SetFormula("1 ++ x");
321 for (Int_t i=0;i<160;++i) {
322 AliTPCclusterMI *c=track->GetClusterPointer(i);
323 if (c && RejectCluster(c)) continue;
324 if (c&&c->GetDetector()==sector) {
325 Double_t x=c->GetX();
326 fy.AddPoint(&x,c->GetY());
327 fz.AddPoint(&x,c->GetZ());
328 xmax=TMath::Max(xmax,x);
329 xmin=TMath::Min(xmin,x);
334 Double_t a[3]={fy.GetParameter(0),
336 type==kQuadratic?fy.GetParameter(2):0.};
337 Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
338 fy.GetCovarianceMatrixElement(1,0),
339 fy.GetCovarianceMatrixElement(1,1),
340 type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0.,
341 type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0.,
342 type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.};
343 for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters;
344 Double_t b[2]={fz.GetParameter(0),
346 Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
347 fz.GetCovarianceMatrixElement(1,0),
348 fz.GetCovarianceMatrixElement(1,1)};
349 for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters;
352 Double_t alpha=track->GetAlpha();
353 Quadratic2Helix(a,ca,b,cb,0.,p,c);
354 fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
355 Quadratic2Helix(a,ca,b,cb,xmin,p,c);
356 fInner=new AliExternalTrackParam(xmin,alpha,p,c);
357 Quadratic2Helix(a,ca,b,cb,xmax,p,c);
358 fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
361 void AliTPCTracklet::Quadratic2Helix(Double_t *a,Double_t *ca,
362 Double_t *b,Double_t *cb,
364 Double_t *p,Double_t *c) {
365 // y(x)=a[0]+a[1]*x+a[2]*x^2
367 // parametrises the corosponding helix at x0
369 // get the polynoms at x0
370 Double_t a0=x0*x0*a[2] + x0*a[1] + a[0];
371 Double_t a1=2.*x0*a[2] + a[1];
373 Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5])));
374 Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5]));
375 Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]);
376 Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]);
377 Double_t ca21=ca[3]+x0*2.*ca[5];
380 Double_t b0=x0*b[1] + b[0];
382 Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]);
383 Double_t cb10=cb[1]+x0*cb[2];
386 // transform to helix parameters
387 Double_t f =1.+a1*a1;
390 Double_t fi12=TMath::Sqrt(fi);
391 Double_t fi32=fi*fi12;
393 Double_t fi52=fi2*fi12;
394 Double_t fi3 =fi2*fi;
395 Double_t fi5 =fi2*fi3;
397 Double_t xyz[3]={0.}; // TODO...
398 Double_t fc=1./(GetBz(xyz)*kB2C);
404 p[4]=2.*a2*fi32*fc; // 1/pt
409 c[3] =ca10*fi32; // snp-y0
411 c[5] =ca11*fi3; // snp-snp
413 c[7] =cb10; // tgl-z0
415 c[9] =cb11; // tgl-tgl
416 c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc; // 1/pt-y0
418 c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp
419 c[13]=0.; // 1/pt-tgl
420 c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5
425 void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) {
427 fy.StoreData(kFALSE);
428 fy.SetFormula("hyp2");
431 for (Int_t i=0;i<160;++i) {
432 AliTPCclusterMI *c=track->GetClusterPointer(i);
433 if (c && RejectCluster(c)) continue;
434 if (c&&c->GetDetector()==sector) {
435 Double_t x=c->GetX();
436 Double_t y=c->GetY();
437 Double_t xy[2]={x,y};
439 Double_t errx=1.,erry=1.;//TODO!
440 Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry);
442 fy.AddPoint(xy,r,err);
443 xmax=TMath::Max(xmax,x);
444 xmin=TMath::Min(xmin,x);
448 Double_t a[3]={fy.GetParameter(0),
451 Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
452 fy.GetCovarianceMatrixElement(1,0),
453 fy.GetCovarianceMatrixElement(1,1),
454 fy.GetCovarianceMatrixElement(2,0),
455 fy.GetCovarianceMatrixElement(2,1),
456 fy.GetCovarianceMatrixElement(2,2)};
459 fz.StoreData(kFALSE);
460 fz.SetFormula("hyp1");
461 Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
465 for (Int_t i=0;i<160;++i) {
466 AliTPCclusterMI *c=track->GetClusterPointer(i);
467 if (c && RejectCluster(c)) continue;
468 if (c&&c->GetDetector()==sector) {
469 Double_t x=c->GetX();
470 Double_t y=c->GetY();
473 phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R));
475 fz.AddPoint(&phi,c->GetZ(),err);
481 Double_t b[2]={fz.GetParameter(0),
483 Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
484 fz.GetCovarianceMatrixElement(1,0),
485 fz.GetCovarianceMatrixElement(1,1)};
489 Double_t alpha=track->GetAlpha();
490 if (Riemann2Helix(a,ca,b,cb,0.,p,c))
491 fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
492 if (Riemann2Helix(a,ca,b,cb,xmin,p,c))
493 fInner=new AliExternalTrackParam(xmin,alpha,p,c);
494 if (Riemann2Helix(a,ca,b,cb,xmax,p,c))
495 fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
498 Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/,
499 Double_t *b,Double_t */*cb*/,
501 Double_t *p,Double_t *c) {
504 Double_t xr0=.5*a[1];
505 Double_t yr0=.5*a[2];
506 Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
508 if (dx*dx>=R*R) return kFALSE;
509 Double_t dy=TMath::Sqrt((R-dx)*(R+dx)); //sign!!
510 if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
513 Double_t tgp=-dx/dy; //TODO: dy!=0
514 Double_t z0=b[0]+TMath::ATan(tgp)*b[1];
515 Double_t xyz[3]={x0,y0,z0};
516 Double_t fc=1./(GetBz(xyz)*kB2C);
520 p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp
522 p[4]=1./R*fc; // 1/pt
536 c[12]=0.; // 1/pt-snp
537 c[13]=0.; // 1/pt-tgl
538 c[14]=0.; // 1/pt-1/pt
543 TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
545 Bool_t storeClusters,
547 Int_t maxTracklets) {
548 // The tracklet factory: It creates several tracklets out of a track. They
549 // are created for sectors that fullfill the constraint of having enough
550 // clusters inside. Futhermore you can specify the maximum amount of
551 // tracklets that are to be created.
552 // The tracklets appear in a sorted fashion, beginning with those having the
555 Int_t sectors[72]={0};
556 for (Int_t i=0;i<160;++i) {
557 AliTPCclusterMI *c=track->GetClusterPointer(i);
558 if (c && RejectCluster(c)) continue;
560 ++sectors[c->GetDetector()];
563 TMath::Sort(72,sectors,indices);
565 if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
566 for (Int_t i=0;i<maxTracklets&§ors[indices[i]]>=minClusters;++i) {
567 tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
572 TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &/*clusters*/,
574 Bool_t /*storeClusters*/,
575 Int_t /*minClusters*/,
576 Int_t /*maxTracklets*/) {
583 Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
584 const AliTPCTracklet &t2,
585 AliExternalTrackParam *&t1m,
586 AliExternalTrackParam *&t2m) {
587 // This function propagates two Tracklets to a common x-coordinate. This
588 // x is dermined as the one that is in the middle of the two tracklets (they
589 // are assumed to live on two distinct x-intervalls).
590 // The inner parametrisation of the outer Tracklet and the outer
591 // parametrisation of the inner Tracklet are used and propagated to this
592 // common x. This result is saved not inside the Tracklets but two new
593 // ExternalTrackParams are created (that means you might want to delete
595 // In the case that the alpha angles of the Tracklets differ both angles
596 // are tried out for this propagation.
597 // In case of any failure kFALSE is returned, no AliExternalTrackParam
598 // is created und the pointers are set to 0.
600 if (t1.GetInner() && t1.GetOuter() &&
601 t2.GetInner() && t2.GetOuter()) {
602 if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) {
603 t1m=new AliExternalTrackParam(*t1.GetOuter());
604 t2m=new AliExternalTrackParam(*t2.GetInner());
607 t1m=new AliExternalTrackParam(*t1.GetInner());
608 t2m=new AliExternalTrackParam(*t2.GetOuter());
610 Double_t mx=.5*(t1m->GetX()+t2m->GetX());
615 Double_t b1[3]; AliTracker::GetBxByBz(xyz,b1);
618 Double_t b2[3]; AliTracker::GetBxByBz(xyz,b2);
619 if (t1m->Rotate(t2m->GetAlpha())
620 //&& t1m->PropagateTo(mx,b1)
621 //&& t2m->PropagateTo(mx,b2));
622 && t1m->PropagateToBxByBz(mx,b1)
623 && t2m->PropagateToBxByBz(mx,b2));
625 if (t2m->Rotate(t1m->GetAlpha())
626 //&& t1m->PropagateTo(mx,b1)
627 //&& t2m->PropagateTo(mx,b2));
628 && t1m->PropagateToBxByBz(mx,b1)
629 && t2m->PropagateToBxByBz(mx,b2));
642 double AliTPCTracklet::GetBz(Double_t *xyz)
644 return AliTracker::GetBz(xyz);
647 void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c,
649 // This function generates a n-dimensional random variable x with mean
650 // p and covariance c.
651 // That is done using the cholesky decomposition of the covariance matrix,
652 // Begin_Latex C=U^{t} U End_Latex, with Begin_Latex U End_Latex being an
653 // upper triangular matrix. Given a vector v of iid gausian random variables
654 // with variance 1 one obtains the asked result as: Begin_Latex x=U^t v
656 // c is expected to be in a lower triangular format:
661 static TRandom3 random;
662 Double_t *c2= new Double_t[ndim*ndim];
664 for (Int_t i=0;i<ndim;++i)
665 for (Int_t j=0;j<=i;++j)
666 c2[i*ndim+j]=c2[j*ndim+i]=c[k++];
667 TMatrixDSym cm(ndim,c2);
669 TDecompChol chol(cm);
671 const TVectorD pv(ndim);
672 const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p));
675 for (Int_t i=0;i<ndim;++i)
677 TMatrixD L=chol.GetU();
682 TEllipse AliTPCTracklet::ErrorEllipse(Double_t x,Double_t y,
683 Double_t sx,Double_t sy,Double_t sxy) {
685 r_{1,2}=1/2 (a+c#pm#sqrt{(a-c)^{2}+(2b)^{2}})
687 Double_t det1=1./(sx*sy-sxy*sxy);
689 Double_t b=-sxy*det1;
692 Double_t s=TMath::Sqrt(d*d+4.*b*b);
693 Double_t r1=TMath::Sqrt(.5*(a+c-s));
694 Double_t r2=TMath::Sqrt(.5*(a+c+s));
695 Double_t alpha=.5*TMath::ATan2(2.*b,d);
696 return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg());
699 void AliTPCTracklet::Test(const char* filename) {
702 AliTPCTracklet::Test("");
703 TFile f("AliTPCTrackletDebug.root");
704 TTree *t=f.Get("AliTPCTrackletDebug");
706 TEllipse e=AliTPCTracklet::ErrorEllipse(0.,0.,4.,1.,1.8);
709 TTreeSRedirector ds(filename);
716 for (Int_t i=0;i<10000;++i) {
719 ds<<"AliTPCTrackletDebug"
732 Double_t param[5]={0.};
733 Double_t covar[15]={1.,
738 AliExternalTrackParam track(x,alpha,param,covar);
742 for (Int_t i=0;i<points.GetNPoints();++i) {
745 Double_t param[5]={0.};
746 Double_t covar[15]={1.,
751 AliExternalTrackParam track(x,alpha,param,covar);
752 for (x=90.;x<250.;x+=1.) {
753 track.PropagateTo(x,b);
761 Bool_t AliTPCTracklet::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
763 // check the acceptance of cluster
764 // Cut on edge effects
766 Bool_t isReject = kFALSE;
767 Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
768 Float_t dist = edgeY - TMath::Abs(cl->GetY());
769 if (param) dist = edgeY - TMath::Abs(param->GetY());
770 if (dist<fgEdgeCutY) isReject=kTRUE;
771 if (cl->GetType()<0) isReject=kTRUE;