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
  41 AliTPCTracklet::AliTPCTracklet()
  42   : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0),
  43     fInner(0),fPrimary(0) {
  44   ////
  45   // The default constructor. It is intended to be used for I/O only.
  46   ////
  47 }
  48
  49 AliTPCTracklet::AliTPCTracklet(const AliTPCseed *track,Int_t sector,
  50                                TrackType type,Bool_t storeClusters)
  51   : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0),
  52     fInner(0),fPrimary(0) {
  53   ////
  54   // Contructor for a tracklet out of a track. Only clusters within a given
  55   // sector are used.
  56   ///
  57
  58   //TODO: only kalman works
  59
  60   for (Int_t i=0;i<160;++i) {
  61     AliTPCclusterMI *c=track->GetClusterPointer(i);
  62     if (c&&c->GetDetector()==sector)
  63       ++fNClusters;
  64   }
  65
  66   if (storeClusters) {
  67     fClusters=new AliTPCclusterMI[fNClusters];
  68     for (Int_t i=0;i<160;++i) {
  69       AliTPCclusterMI *c=track->GetClusterPointer(i);
  70       if (c&&c->GetDetector()==sector)
  71         fClusters[fNStoredClusters]=c;
  72       ++fNStoredClusters;
  73     }
  74   }
  75
  76   switch (type) {
  77   case kKalman:
  78     FitKalman(track,sector);
  79     break;
  80   case kLinear:
  81   case kQuadratic:
  82     FitLinear(track,sector,type);
  83     break;
  84   case kRiemann:
  85     FitRiemann(track,sector);
  86     break;
  87   }
  88
  89 }
  90
  91 AliTPCTracklet::AliTPCTracklet(const TObjArray &clusters,Int_t sector,
  92                                TrackType type,Bool_t storeClusters) {
  93   //TODO: write it!
  94 }
  95
  96 AliTPCTracklet::AliTPCTracklet(const AliTPCTracklet &t)
  97   : fNClusters(t.fNClusters),fNStoredClusters(t.fNStoredClusters),fClusters(0),
  98     fSector(t.fSector),fOuter(0),fInner(0),
  99     fPrimary(0) {
 100   ////
 101   // The copy constructor. You can copy tracklets!
 102   ////
 103
 104   if (t.fClusters) {
 105     fClusters=new AliTPCclusterMI[t.fNStoredClusters];
 106     for (int i=0;i<t.fNStoredClusters;++i)
 107       fClusters[i]=t.fClusters[i];
 108   }
 109   if (t.fOuter)
 110     fOuter=new AliExternalTrackParam(*t.fOuter);
 111   if (t.fInner)
 112     fInner=new AliExternalTrackParam(*t.fInner);
 113   if (t.fPrimary)
 114     fPrimary=new AliExternalTrackParam(*t.fPrimary);
 115 }
 116
 117 AliTPCTracklet& AliTPCTracklet::operator=(const AliTPCTracklet &t) {
 118   ////
 119   // The assignment constructor. You can assign tracklets!
 120   ////
 121   if (this!=&t) {
 122     fNClusters=t.fNClusters;
 123     fNStoredClusters=fNStoredClusters;
 124     delete fClusters;
 125     if (t.fClusters) {
 126       fClusters=new AliTPCclusterMI[t.fNStoredClusters];
 127       for (int i=0;i<t.fNStoredClusters;++i)
 128         fClusters[i]=t.fClusters[i];
 129     }
 130     else
 131       fClusters=0;
 132     fSector=t.fSector;
 133     if (t.fOuter) {
 134       if (fOuter)
 135         *fOuter=*t.fOuter;
 136       else
 137         fOuter=new AliExternalTrackParam(*t.fOuter);
 138     }
 139     else {
 140       delete fOuter;
 141       fOuter=0;
 142     }
 143
 144     if (t.fInner) {
 145       if (fInner)
 146         *fInner=*t.fInner;
 147       else
 148         fInner=new AliExternalTrackParam(*t.fInner);
 149     }
 150     else {
 151       delete fInner;
 152       fInner=0;
 153     }
 154
 155     if (t.fPrimary) {
 156       if (fPrimary)
 157         *fPrimary=*t.fPrimary;
 158       else
 159         fPrimary=new AliExternalTrackParam(*t.fPrimary);
 160     }
 161     else {
 162       delete fPrimary;
 163       fPrimary=0;
 164     }
 165   }
 166   return *this;
 167 }
 168
 169 AliTPCTracklet::~AliTPCTracklet() {
 170   //
 171   // The destructor. Yes, you can even destruct tracklets.
 172   //
 173   delete fClusters;
 174   delete fOuter;
 175   delete fInner;
 176   delete fPrimary;
 177 }
 178
 179 void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
 180   AliTPCseed *t=new AliTPCseed(*track);
 181   if (!t->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-t->GetAlpha())) {
 182     delete t;
 183     return;
 184   }
 185   // fit from inner to outer row
 186   AliTPCseed *outerSeed=new AliTPCseed(*t);
 187   Int_t n=0;
 188   for (Int_t i=0;i<160;++i) {
 189     AliTPCclusterMI *c=t->GetClusterPointer(i);
 190     if (c&&c->GetDetector()==sector) {
 191       if (n==1) {
 192         outerSeed->ResetCovariance(100.);
 193       }
 194       ++n;
 195       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 196       Double_t cov[3]={0.1,0.,0.1}; //TODO: correct error parametrisation
 197       if (!outerSeed->PropagateTo(r[0]) ||
 198           !static_cast<AliExternalTrackParam*>(outerSeed)->Update(&r[1],cov)) {
 199         delete outerSeed;
 200         outerSeed=0;
 201         break;
 202       }
 203     }
 204   }
 205   if (outerSeed)
 206     fOuter=new AliExternalTrackParam(*outerSeed);
 207   delete outerSeed;
 208   // fit from outer to inner rows
 209   AliTPCseed *innerSeed=new AliTPCseed(*t);
 210   n=0;
 211   for (Int_t i=159;i>=0;--i) {
 212     AliTPCclusterMI *c=t->GetClusterPointer(i);
 213     if (c&&c->GetDetector()==sector) {
 214       if (n==1) {
 215         innerSeed->ResetCovariance(100.);
 216       }
 217       ++n;
 218       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 219       Double_t cov[3]={0.1,0.,0.1};
 220       if (!innerSeed->PropagateTo(r[0]) ||
 221           !static_cast<AliExternalTrackParam*>(innerSeed)->Update(&r[1],cov)) {
 222         delete innerSeed;
 223         innerSeed=0;
 224         break;
 225       }
 226     }
 227   }
 228   if (innerSeed)
 229     fInner=new AliExternalTrackParam(*innerSeed);
 230   // propagate to the primary vertex
 231   if (innerSeed) {
 232     AliTPCseed *primarySeed=new AliTPCseed(*innerSeed);
 233     Double_t pos[]={0.,0.,0.};
 234     Double_t sigma[]={.1,.1,.1}; //TODO: is this correct?
 235     AliESDVertex vertex(pos,sigma);
 236     if (primarySeed->PropagateToVertex(&vertex))
 237       fPrimary=new AliExternalTrackParam(*primarySeed);
 238     delete primarySeed;
 239     // for better comparison one does not want to have alpha changed...
 240     if (!fPrimary->Rotate(fInner->GetAlpha())) {
 241       delete fPrimary;
 242       fPrimary=0;
 243     }
 244   }
 245   delete innerSeed;
 246
 247   delete t;
 248 }
 249
 250 void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
 251                                TrackType type) {
 252   TLinearFitter fy(1);
 253   TLinearFitter fz(1);
 254   fy.StoreData(kFALSE);
 255   fz.StoreData(kFALSE);
 256   switch (type) {
 257   case kLinear:
 258     fy.SetFormula("1 ++ x");
 259     fz.SetFormula("1 ++ x");
 260     break;
 261   case kQuadratic:
 262     fy.SetFormula("1 ++ x ++ x*x");
 263     fz.SetFormula("1 ++ x");
 264     break;
 265   }
 266   Double_t xmax=-1.;
 267   Double_t xmin=1000.;
 268   for (Int_t i=0;i<160;++i) {
 269     AliTPCclusterMI *c=track->GetClusterPointer(i);
 270     if (c&&c->GetDetector()==sector) {
 271       Double_t x=c->GetX();
 272       fy.AddPoint(&x,c->GetY());
 273       fz.AddPoint(&x,c->GetZ());
 274       xmax=TMath::Max(xmax,x);
 275       xmin=TMath::Min(xmin,x);
 276     }
 277   }
 278   fy.Eval();
 279   fz.Eval();
 280   Double_t a[3]={fy.GetParameter(0),
 281                  fy.GetParameter(1),
 282                  type==kQuadratic?fy.GetParameter(2):0.};
 283   Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
 284                   fy.GetCovarianceMatrixElement(1,0),
 285                   fy.GetCovarianceMatrixElement(1,1),
 286                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0.,
 287                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0.,
 288                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.};
 289   for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters;
 290   Double_t b[2]={fz.GetParameter(0),
 291                  fz.GetParameter(1)};
 292   Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
 293                   fz.GetCovarianceMatrixElement(1,0),
 294                   fz.GetCovarianceMatrixElement(1,1)};
 295   for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters;
 296   Double_t p[5];
 297   Double_t c[15];
 298   Double_t alpha=track->GetAlpha();
 299   Quadratic2Helix(a,ca,b,cb,0.,p,c);
 300   fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
 301   Quadratic2Helix(a,ca,b,cb,xmin,p,c);
 302   fInner=new AliExternalTrackParam(xmin,alpha,p,c);
 303   Quadratic2Helix(a,ca,b,cb,xmax,p,c);
 304   fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
 305 }
 306
 307 void AliTPCTracklet::Quadratic2Helix(Double_t *a,Double_t *ca,
 308                                      Double_t *b,Double_t *cb,
 309                                      Double_t x0,
 310                                      Double_t *p,Double_t *c) {
 311   // y(x)=a[0]+a[1]*x+a[2]*x^2
 312   // z(x)=b[0]+b[1]*x
 313   // parametrises the corosponding helix at x0
 314
 315   // get the polynoms at x0
 316   Double_t a0=x0*x0*a[2] + x0*a[1] + a[0];
 317   Double_t a1=2.*x0*a[2] +     a[1];
 318   Double_t a2=      a[2];
 319   Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5])));
 320   Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5]));
 321   Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]);
 322   Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]);
 323   Double_t ca21=ca[3]+x0*2.*ca[5];
 324   Double_t ca22=ca[5];
 325
 326   Double_t b0=x0*b[1] + b[0];
 327   Double_t b1=   b[1];
 328   Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]);
 329   Double_t cb10=cb[1]+x0*cb[2];
 330   Double_t cb11=cb[2];
 331
 332   // transform to helix parameters
 333   Double_t f   =1.+a1*a1;
 334   Double_t f2  =f*f;
 335   Double_t fi  =1./f;
 336   Double_t fi12=TMath::Sqrt(fi);
 337   Double_t fi32=fi*fi12;
 338   Double_t fi2 =fi*fi;
 339   Double_t fi52=fi2*fi12;
 340   Double_t fi3 =fi2*fi;
 341   Double_t fi5 =fi2*fi3;
 342
 343   Double_t xyz[3]={0.}; // TODO...
 344   Double_t fc=1./(GetBz(xyz)*kB2C);
 345
 346   p[0]=a0;            // y0
 347   p[1]=b0;            // z0
 348   p[2]=a1*fi12;       // snp
 349   p[3]=b1;            // tgl
 350   p[4]=2.*a2*fi32*fc; // 1/pt
 351
 352   c[0] =ca00;      //  y0-y0
 353   c[1] =0.;        //  z0-y0
 354   c[2] =cb00;      //  z0-z0
 355   c[3] =ca10*fi32; // snp-y0
 356   c[4] =0.;        // snp-z0
 357   c[5] =ca11*fi3;  // snp-snp
 358   c[6] =0.;        // tgl-y0
 359   c[7] =cb10;      // tgl-z0
 360   c[8] =0.;        // tgl-snp
 361   c[9] =cb11;      // tgl-tgl
 362   c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc;  // 1/pt-y0
 363   c[11]=0.;                                 // 1/pt-z0
 364   c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp
 365   c[13]=0.;                                 // 1/pt-tgl
 366   c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5
 367     *fc*fc;        // 1/pt-1/pt
 368 }
 369
 370
 371 void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) {
 372   TLinearFitter fy(2);
 373   fy.StoreData(kFALSE);
 374   fy.SetFormula("hyp2");
 375   Double_t xmax=-1.;
 376   Double_t xmin=1000.;
 377   for (Int_t i=0;i<160;++i) {
 378     AliTPCclusterMI *c=track->GetClusterPointer(i);
 379     if (c&&c->GetDetector()==sector) {
 380       Double_t x=c->GetX();
 381       Double_t y=c->GetY();
 382       Double_t xy[2]={x,y};
 383       Double_t r=x*x+y*y;
 384       Double_t errx=1.,erry=1.;//TODO!
 385       Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry);
 386       err=1.;
 387       fy.AddPoint(xy,r,err);
 388       xmax=TMath::Max(xmax,x);
 389       xmin=TMath::Min(xmin,x);
 390     }
 391   }
 392   fy.Eval();
 393   Double_t a[3]={fy.GetParameter(0),
 394                  fy.GetParameter(1),
 395                  fy.GetParameter(2)};
 396   Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
 397                   fy.GetCovarianceMatrixElement(1,0),
 398                   fy.GetCovarianceMatrixElement(1,1),
 399                   fy.GetCovarianceMatrixElement(2,0),
 400                   fy.GetCovarianceMatrixElement(2,1),
 401                   fy.GetCovarianceMatrixElement(2,2)};
 402
 403   TLinearFitter fz(1);
 404   fz.StoreData(kFALSE);
 405   fz.SetFormula("hyp1");
 406   Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
 407   Double_t oldx=0.;
 408   Double_t oldy=R;
 409   Double_t phi=0.;
 410   for (Int_t i=0;i<160;++i) {
 411     AliTPCclusterMI *c=track->GetClusterPointer(i);
 412     if (c&&c->GetDetector()==sector) {
 413       Double_t x=c->GetX();
 414       Double_t y=c->GetY();
 415       Double_t dx=x-oldx;
 416       Double_t dy=y-oldy;
 417       phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R));
 418       Double_t err=1.;
 419       fz.AddPoint(&phi,c->GetZ(),err);
 420       oldx=x;
 421       oldy=y;
 422     }
 423   }
 424   fz.Eval();
 425   Double_t b[2]={fz.GetParameter(0),
 426                  fz.GetParameter(1)};
 427   Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
 428                   fz.GetCovarianceMatrixElement(1,0),
 429                   fz.GetCovarianceMatrixElement(1,1)};
 430
 431   Double_t p[5];
 432   Double_t c[15];
 433   Double_t alpha=track->GetAlpha();
 434   if (Riemann2Helix(a,ca,b,cb,0.,p,c))
 435     fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
 436   if (Riemann2Helix(a,ca,b,cb,xmin,p,c))
 437     fInner=new AliExternalTrackParam(xmin,alpha,p,c);
 438   if (Riemann2Helix(a,ca,b,cb,xmax,p,c))
 439     fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
 440 }
 441
 442 Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t *ca,
 443                                      Double_t *b,Double_t *cb,
 444                                      Double_t x0,
 445                                      Double_t *p,Double_t *c) {
 446   //TODO: signs!
 447
 448   Double_t xr0=.5*a[1];
 449   Double_t yr0=.5*a[2];
 450   Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
 451   Double_t dx=x0-xr0;
 452   if (dx*dx>=R*R) return kFALSE;
 453   Double_t dy=TMath::Sqrt(R*R-dx*dx); //sign!!
 454   if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
 455     dy=-dy;
 456   Double_t y0=yr0+dy;
 457   Double_t tgp=-dx/dy; //TODO: dy!=0
 458   Double_t z0=b[0]+TMath::ATan(tgp)*b[1];
 459   Double_t xyz[3]={x0,y0,z0};
 460   Double_t fc=1./(GetBz(xyz)*kB2C);
 461   fc=1;
 462   p[0]=y0;  // y0
 463   p[1]=z0; // z0
 464   p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp
 465   p[3]=b[1];       // tgl
 466   p[4]=1./R*fc;    // 1/pt
 467
 468   c[0] =0.;      //  y0-y0
 469   c[1] =0.;        //  z0-y0
 470   c[2] =0.;      //  z0-z0
 471   c[3] =0.; // snp-y0
 472   c[4] =0.;        // snp-z0
 473   c[5] =0.;  // snp-snp
 474   c[6] =0.;        // tgl-y0
 475   c[7] =0.;      // tgl-z0
 476   c[8] =0.;        // tgl-snp
 477   c[9] =0.;      // tgl-tgl
 478   c[10]=0.;  // 1/pt-y0
 479   c[11]=0.;                                 // 1/pt-z0
 480   c[12]=0.; // 1/pt-snp
 481   c[13]=0.;                                 // 1/pt-tgl
 482   c[14]=0.;        // 1/pt-1/pt
 483
 484   return kTRUE;
 485 }
 486
 487 TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
 488                                           TrackType type,
 489                                           Bool_t storeClusters,
 490                                           Int_t minClusters,
 491                                           Int_t maxTracklets) {
 492 // The tracklet factory: It creates several tracklets out of a track. They
 493 // are created for sectors that fullfill the constraint of having enough
 494 // clusters inside. Futhermore you can specify the maximum amount of
 495 // tracklets that are to be created.
 496 // The tracklets appear in a sorted fashion, beginning with those having the
 497 // most clusters.
 498
 499   Int_t sectors[72]={0};
 500   for (Int_t i=0;i<160;++i) {
 501     AliTPCclusterMI *c=track->GetClusterPointer(i);
 502     if (c)
 503       ++sectors[c->GetDetector()];
 504   }
 505   Int_t indices[72];
 506   TMath::Sort(72,sectors,indices);
 507   TObjArray tracklets;
 508   if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
 509   for (Int_t i=0;i<maxTracklets&&sectors[indices[i]]>=minClusters;++i) {
 510     tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
 511   }
 512   return tracklets;
 513 }
 514
 515 TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &clusters,
 516                                           TrackType type,
 517                                           Bool_t storeClusters,
 518                                           Int_t minClusters,
 519                                           Int_t maxTracklets) {
 520   TObjArray tracklets;
 521   return tracklets;
 522 }
 523
 524 Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
 525                                         const AliTPCTracklet &t2,
 526                                         AliExternalTrackParam *&t1m,
 527                                         AliExternalTrackParam *&t2m) {
 528   // This function propagates two Tracklets to a common x-coordinate. This
 529   // x is dermined as the one that is in the middle of the two tracklets (they
 530   // are assumed to live on two distinct x-intervalls).
 531   // The inner parametrisation of the outer Tracklet and the outer
 532   // parametrisation of the inner Tracklet are used and propagated to this
 533   // common x. This result is saved not inside the Tracklets but two new
 534   // ExternalTrackParams are created (that means you might want to delete
 535   // them).
 536   // In the case that the alpha angles of the Tracklets differ both angles
 537   // are tried out for this propagation.
 538   // In case of any failure kFALSE is returned, no AliExternalTrackParam
 539   // is created und the pointers are set to 0.
 540
 541   if (t1.GetInner() && t1.GetOuter() &&
 542       t2.GetInner() && t2.GetOuter()) {
 543     if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) {
 544       t1m=new AliExternalTrackParam(*t1.GetOuter());
 545       t2m=new AliExternalTrackParam(*t2.GetInner());
 546     }
 547     else {
 548       t1m=new AliExternalTrackParam(*t1.GetInner());
 549       t2m=new AliExternalTrackParam(*t2.GetOuter());
 550     }
 551     Double_t mx=.5*(t1m->GetX()+t2m->GetX());
 552     Double_t b1,b2;
 553     Double_t xyz[3];
 554     t1m->GetXYZ(xyz);
 555     b1=GetBz(xyz);
 556     t2m->GetXYZ(xyz);
 557     b2=GetBz(xyz);
 558     if (t1m->Rotate(t2m->GetAlpha())
 559         && t1m->PropagateTo(mx,b1)
 560         && t2m->PropagateTo(mx,b2));
 561     else
 562       if (t2m->Rotate(t1m->GetAlpha())
 563           && t1m->PropagateTo(mx,b1)
 564           && t2m->PropagateTo(mx,b2));
 565       else {
 566         delete t1m;
 567         delete t2m;
 568         t1m=t2m=0;
 569       }
 570   }
 571   else {
 572     t1m=t2m=0;
 573   }
 574   return t1m&&t2m;
 575 }
 576
 577 double AliTPCTracklet::GetBz(Double_t *xyz) {
 578   if (AliTracker::UniformField())
 579     return AliTracker::GetBz();
 580   else
 581     return AliTracker::GetBz(xyz);
 582 }
 583
 584 void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c,
 585                               Double_t *x) {
 586   // This function generates a n-dimensional random variable x with mean
 587   // p and covariance c.
 588   // That is done using the cholesky decomposition of the covariance matrix,
 589   // Begin_Latex C=U^{t} U End_Latex, with Begin_Latex U End_Latex being an
 590   // upper triangular matrix. Given a vector v of iid gausian random variables
 591   // with variance 1 one obtains the asked result as: Begin_Latex x=U^t v
 592   // End_Latex.
 593   // c is expected to be in a lower triangular format:
 594   // c[0]
 595   // c[1] c[2]
 596   // c[3] c[4] c[5]
 597   // etc.
 598   static TRandom3 random;
 599   Double_t *c2= new Double_t[ndim*ndim];
 600   Int_t k=0;
 601   for (Int_t i=0;i<ndim;++i)
 602     for (Int_t j=0;j<=i;++j)
 603       c2[i*ndim+j]=c2[j*ndim+i]=c[k++];
 604   TMatrixDSym cm(ndim,c2);
 605   delete[] c2;
 606   TDecompChol chol(cm);
 607   chol.Decompose();
 608   const TVectorD pv(ndim);
 609   const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p));
 610   TVectorD xv(ndim);
 611   xv.Use(ndim,x);
 612   for (Int_t i=0;i<ndim;++i)
 613     xv[i]=random.Gaus();
 614   TMatrixD L=chol.GetU();
 615   L.T();
 616   xv=L*xv+pv;
 617 }
 618
 619 TEllipse AliTPCTracklet::ErrorEllipse(Double_t x,Double_t y,
 620                                       Double_t sx,Double_t sy,Double_t sxy) {
 621   /* Begin_Latex
 622      r_{1,2}=1/2 (a+c#pm#sqrt{(a-c)^{2}+(2b)^{2}})
 623   End_Latex */
 624   Double_t det1=1./(sx*sy-sxy*sxy);
 625   Double_t a=sy*det1;
 626   Double_t b=-sxy*det1;
 627   Double_t c=sx*det1;
 628   Double_t d=c-a;
 629   Double_t s=TMath::Sqrt(d*d+4.*b*b);
 630   Double_t r1=TMath::Sqrt(.5*(a+c-s));
 631   Double_t r2=TMath::Sqrt(.5*(a+c+s));
 632   Double_t alpha=.5*TMath::ATan2(2.*b,d);
 633   return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg());
 634 }
 635
 636 void AliTPCTracklet::Test(const char* filename) {
 637   /*
 638     aliroot
 639     AliTPCTracklet::Test("");
 640     TFile f("AliTPCTrackletDebug.root");
 641     TTree *t=f.Get("AliTPCTrackletDebug");
 642     t->Draw("p0:p4");
 643     TEllipse e=AliTPCTracklet::ErrorEllipse(0.,0.,4.,1.,1.8);
 644     e.Draw();
 645  */
 646   TTreeSRedirector ds("AliTPCTrackletDebug.root");
 647   Double_t p[5]={0.};
 648   Double_t c[15]={4.,
 649                   0.,4.,
 650                   0.,0.,9.,
 651                   0.,0.,0.,16.,
 652                   1.8,0.,0.,0.,1.};
 653   for (Int_t i=0;i<10000;++i) {
 654     Double_t x[5];
 655     RandomND(5,p,c,x);
 656     ds<<"AliTPCTrackletDebug"
 657       <<"p0="<<x[0]
 658       <<"p1="<<x[1]
 659       <<"p2="<<x[2]
 660       <<"p3="<<x[3]
 661       <<"p4="<<x[4]
 662       <<"\n";
 663   }
 664
 665   /*
 666   Double_t b;
 667   Double_t x=0.;
 668   Double_t alpha=0.;
 669   Double_t param[5]={0.};
 670   Double_t covar[15]={1.,
 671                       0.,1.,
 672                       0.,0.,1.,
 673                       0.,0.,0.,1.,
 674                       0.,0.,0.,0.,1.};
 675   AliExternalTrackParam track(x,alpha,param,covar);
 676
 677
 678
 679   for (Int_t i=0;i<points.GetNPoints();++i) {
 680     Double_t x=0.;
 681     Double_t alpha=0.;
 682     Double_t param[5]={0.};
 683     Double_t covar[15]={1.,
 684                         0.,1.,
 685                         0.,0.,1.,
 686                         0.,0.,0.,1.,
 687                         0.,0.,0.,0.,1.};
 688     AliExternalTrackParam track(x,alpha,param,covar);
 689     for (x=90.;x<250.;x+=1.) {
 690       track.PropagateTo(x,b);
 691       AliTPCclusterMI c();
 692     }
 693   }
 694   */
 695 }