TPC/AliTPCTracklet.cxx

   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
  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"
  31 #include "TRandom3.h"
  32 #include "TDecompChol.h"
  33
  34 #include <iostream>
  35 using namespace std;
  36
  37 ClassImp(AliTPCTracklet)
  38
  39 const Double_t AliTPCTracklet::kB2C=0.299792458e-3;
  40 Float_t  AliTPCTracklet::fgEdgeCutY=3;
  41 Float_t  AliTPCTracklet::fgEdgeCutX=0;
  42
  43 AliTPCTracklet::AliTPCTracklet()
  44   : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0),
  45     fInner(0),fPrimary(0) {
  46   ////
  47   // The default constructor. It is intended to be used for I/O only.
  48   ////
  49 }
  50
  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) {
  55   ////
  56   // Contructor for a tracklet out of a track. Only clusters within a given
  57   // sector are used.
  58   ///
  59
  60   //TODO: only kalman works
  61
  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)
  66       ++fNClusters;
  67   }
  68
  69   if (storeClusters) {
  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;
  76       ++fNStoredClusters;
  77     }
  78   }
  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;
  91   }
  92
  93 }
  94
  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) {
  99   //TODO: write it!
 100 }
 101
 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),
 105     fPrimary(0) {
 106   ////
 107   // The copy constructor. You can copy tracklets!
 108   ////
 109
 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   }
 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   ////
 127   if (this!=&t) {
 128     fNClusters=t.fNClusters;
 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;
 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
 174 AliTPCTracklet::~AliTPCTracklet() {
 175   //
 176   // The destructor. Yes, you can even destruct tracklets.
 177   //
 178   delete fClusters;
 179   delete fOuter;
 180   delete fInner;
 181   delete fPrimary;
 182 }
 183
 184
 185
 186
 187
 188 void 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;
 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
 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;
 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
 235     AliTPCseed::GetError(c, track,cov[0],cov[2]);
 236     cov[0]*=cov[0];
 237     cov[2]*=cov[2];
 238     if (!track->PropagateTo(r[0])) {
 239       isOK=kFALSE;
 240       break;
 241     }
 242     if (RejectCluster(c,track)) continue;
 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;
 253     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 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];
 258     if (!track->PropagateTo(r[0])) {
 259       isOK=kFALSE;
 260       break;
 261     }
 262     if (RejectCluster(c,track)) continue;
 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;
 273     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 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];
 278     if (!track->PropagateTo(r[0])) {
 279       isOK=kFALSE;
 280       break;
 281     }
 282     if (RejectCluster(c,track)) continue;
 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   //
 293   delete track;
 294 }
 295
 296
 297
 298
 299 void 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;
 314   case kKalman:
 315   case kRiemann:
 316     break;
 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);
 322     if (c && RejectCluster(c)) continue;
 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
 360 void 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
 424 void 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);
 432     if (c && RejectCluster(c)) continue;
 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);
 466     if (c && RejectCluster(c)) continue;
 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
 497 Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/,
 498                                      Double_t *b,Double_t */*cb*/,
 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;
 508   Double_t dy=TMath::Sqrt((R-dx)*(R+dx)); //sign!!
 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
 542 TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
 543                                           TrackType type,
 544                                           Bool_t storeClusters,
 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) {
 556     AliTPCclusterMI *c=track->GetClusterPointer(i);
 557     if (c && RejectCluster(c)) continue;
 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) {
 566     tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
 567   }
 568   return tracklets;
 569 }
 570
 571 TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &/*clusters*/,
 572                                           TrackType /*type*/,
 573                                           Bool_t /*storeClusters*/,
 574                                           Int_t /*minClusters*/,
 575                                           Int_t /*maxTracklets*/) {
 576   // TODO!
 577
 578   TObjArray tracklets;
 579   return tracklets;
 580 }
 581
 582 Bool_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());
 610     //Double_t b1,b2;
 611     Double_t xyz[3];
 612     t1m->GetXYZ(xyz);
 613     //b1=GetBz(xyz);
 614     Double_t b1[3]; AliTracker::GetBxByBz(xyz,b1);
 615     t2m->GetXYZ(xyz);
 616     //b2=GetBz(xyz);
 617     Double_t b2[3]; AliTracker::GetBxByBz(xyz,b2);
 618     if (t1m->Rotate(t2m->GetAlpha())
 619         //&& t1m->PropagateTo(mx,b1)
 620         //&& t2m->PropagateTo(mx,b2));
 621         && t1m->PropagateToBxByBz(mx,b1)
 622         && t2m->PropagateToBxByBz(mx,b2));
 623     else
 624       if (t2m->Rotate(t1m->GetAlpha())
 625           //&& t1m->PropagateTo(mx,b1)
 626           //&& t2m->PropagateTo(mx,b2));
 627           && t1m->PropagateToBxByBz(mx,b1)
 628           && t2m->PropagateToBxByBz(mx,b2));
 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
 641 double AliTPCTracklet::GetBz(Double_t *xyz)
 642 {
 643   return AliTracker::GetBz(xyz);
 644 }
 645
 646 void 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
 681 TEllipse 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
 698 void 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  */
 708   TTreeSRedirector ds(filename);
 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 }
 758
 759
 760 Bool_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 }