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