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