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     fNStoredClusters=fNStoredClusters;
 130     delete fClusters;
 131     if (t.fClusters) {
 132       fClusters=new AliTPCclusterMI[t.fNStoredClusters];
 133       for (int i=0;i<t.fNStoredClusters;++i)
 134         fClusters[i]=t.fClusters[i];
 135     }
 136     else
 137       fClusters=0;
 138     fSector=t.fSector;
 139     if (t.fOuter) {
 140       if (fOuter)
 141         *fOuter=*t.fOuter;
 142       else
 143         fOuter=new AliExternalTrackParam(*t.fOuter);
 144     }
 145     else {
 146       delete fOuter;
 147       fOuter=0;
 148     }
 149
 150     if (t.fInner) {
 151       if (fInner)
 152         *fInner=*t.fInner;
 153       else
 154         fInner=new AliExternalTrackParam(*t.fInner);
 155     }
 156     else {
 157       delete fInner;
 158       fInner=0;
 159     }
 160
 161     if (t.fPrimary) {
 162       if (fPrimary)
 163         *fPrimary=*t.fPrimary;
 164       else
 165         fPrimary=new AliExternalTrackParam(*t.fPrimary);
 166     }
 167     else {
 168       delete fPrimary;
 169       fPrimary=0;
 170     }
 171   }
 172   return *this;
 173 }
 174
 175 AliTPCTracklet::~AliTPCTracklet() {
 176   //
 177   // The destructor. Yes, you can even destruct tracklets.
 178   //
 179   delete fClusters;
 180   delete fOuter;
 181   delete fInner;
 182   delete fPrimary;
 183 }
 184
 185
 186
 187 /*
 188 void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
 189   //
 190   // Fit using Kalman filter
 191   //
 192   AliTPCseed *t=new AliTPCseed(*track);
 193   if (!t->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-t->GetAlpha())) {
 194     delete t;
 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]=10.*10.;
 201   covar[2]=10.*10.;
 202   covar[5]=10.*10./(64.*64.);
 203   covar[9]=10.*10./(64.*64.);
 204   covar[14]=1*1;
 205
 206   //
 207   AliTPCseed *outerSeed=new AliTPCseed(*t);
 208   Int_t n=0;
 209   for (Int_t i=0;i<160;++i) {
 210
 211     AliTPCclusterMI *c=t->GetClusterPointer(i);
 212     if (c && RejectCluster(c,outerSeed)) continue;
 213     if (c&&c->GetDetector()==sector) {
 214       if (n==1) {
 215         outerSeed->ResetCovariance(100.);
 216         outerSeed->AddCovariance(covar);
 217       }
 218       ++n;
 219       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 220       Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 221       if (!outerSeed->PropagateTo(r[0]) ||
 222           !static_cast<AliExternalTrackParam*>(outerSeed)->Update(&r[1],cov)) {
 223         delete outerSeed;
 224         outerSeed=0;
 225         break;
 226       }
 227     }
 228   }
 229   if (outerSeed)
 230     fOuter=new AliExternalTrackParam(*outerSeed);
 231   // fit from outer to inner rows
 232   //  AliTPCseed *innerSeed=new AliTPCseed(*t);
 233   AliTPCseed *innerSeed=0;
 234   if (fOuter) innerSeed=new AliTPCseed(*outerSeed);
 235   if (!innerSeed) innerSeed=new AliTPCseed(*t);
 236   delete outerSeed;
 237
 238   n=0;
 239   for (Int_t i=159;i>=0;--i) {
 240     AliTPCclusterMI *c=t->GetClusterPointer(i);
 241     if (c && RejectCluster(c, innerSeed)) continue;
 242     if (c&&c->GetDetector()==sector) {
 243       if (n==1) {
 244         innerSeed->ResetCovariance(100.);
 245         innerSeed->AddCovariance(covar);
 246       }
 247       ++n;
 248       Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 249       Double_t cov[3]={0.01,0.,0.01};
 250       if (!innerSeed->PropagateTo(r[0]) ||
 251           !static_cast<AliExternalTrackParam*>(innerSeed)->Update(&r[1],cov)) {
 252         delete innerSeed;
 253         innerSeed=0;
 254         break;
 255       }
 256     }
 257   }
 258   if (innerSeed)
 259     fInner=new AliExternalTrackParam(*innerSeed);
 260   // propagate to the primary vertex
 261   if (innerSeed) {
 262     AliTPCseed *primarySeed=new AliTPCseed(*innerSeed);
 263     Double_t pos[]={0.,0.,0.};
 264     Double_t sigma[]={.1,.1,.1}; //TODO: is this correct?
 265     AliESDVertex vertex(pos,sigma);
 266     if (primarySeed->PropagateToVertex(&vertex))
 267       fPrimary=new AliExternalTrackParam(*primarySeed);
 268     delete primarySeed;
 269     // for better comparison one does not want to have alpha changed...
 270     if (fPrimary) if (!fPrimary->Rotate(fInner->GetAlpha())) {
 271       delete fPrimary;
 272       fPrimary=0;
 273     }
 274   }
 275   delete innerSeed;
 276
 277   delete t;
 278 }
 279 */
 280
 281
 282 void AliTPCTracklet::FitKalman(const AliTPCseed *seed,Int_t sector) {
 283   //
 284   // Fit using Kalman filter
 285   //
 286   AliTPCseed *track=new AliTPCseed(*seed);
 287   if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
 288     delete track;
 289     return;
 290   }
 291   // fit from inner to outer row
 292   Double_t covar[15];
 293   for (Int_t i=0;i<15;i++) covar[i]=0;
 294   covar[0]=10.*10.;
 295   covar[2]=10.*10.;
 296   covar[5]=10.*10./(64.*64.);
 297   covar[9]=10.*10./(64.*64.);
 298   covar[14]=1*1;
 299   Float_t xmin=1000, xmax=-10000;
 300   Int_t imin=158, imax=0;
 301   for (Int_t i=0;i<160;i++) {
 302     AliTPCclusterMI *c=track->GetClusterPointer(i);
 303     if (!c) continue;
 304     if (c->GetDetector()!=sector)  continue;
 305     if (c->GetX()<xmin) xmin=c->GetX();
 306     if (c->GetX()>xmax) xmax=c->GetX();
 307     if (i<imin) imin=i;
 308     if (i>imax) imax=i;
 309   }
 310   if(imax-imin<10) {
 311     delete track;
 312     return;
 313   }
 314
 315   for (Float_t x=track->GetX(); x<xmin; x++) track->PropagateTo(x);
 316   track->AddCovariance(covar);
 317   //
 318   AliExternalTrackParam paramIn;
 319   AliExternalTrackParam paramOut;
 320   Bool_t isOK=kTRUE;
 321   //
 322   //
 323   //
 324   for (Int_t i=imin; i<=imax; i++){
 325     AliTPCclusterMI *c=track->GetClusterPointer(i);
 326     if (!c) continue;
 327     if (RejectCluster(c,track)) continue;
 328     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 329     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 330     if (!track->PropagateTo(r[0])) {
 331       isOK=kFALSE;
 332       break;
 333     }
 334     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
 335   }
 336   if (!isOK) { delete track; return;}
 337   track->AddCovariance(covar);
 338   //
 339   //
 340   //
 341   for (Int_t i=imax; i>=imin; i--){
 342     AliTPCclusterMI *c=track->GetClusterPointer(i);
 343     if (!c) continue;
 344     if (RejectCluster(c,track)) continue;
 345     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 346     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 347     if (!track->PropagateTo(r[0])) {
 348       isOK=kFALSE;
 349       break;
 350     }
 351     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
 352   }
 353   if (!isOK) { delete track; return;}
 354   paramIn = *track;
 355   track->AddCovariance(covar);
 356   //
 357   //
 358   for (Int_t i=imin; i<=imax; i++){
 359     AliTPCclusterMI *c=track->GetClusterPointer(i);
 360     if (!c) continue;
 361     if (RejectCluster(c,track)) continue;
 362     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
 363     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 364     if (!track->PropagateTo(r[0])) {
 365       isOK=kFALSE;
 366       break;
 367     }
 368     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
 369   }
 370   if (!isOK) { delete track; return;}
 371   paramOut=*track;
 372   //
 373   //
 374   //
 375   fOuter=new AliExternalTrackParam(paramOut);
 376   fInner=new AliExternalTrackParam(paramIn);
 377   //
 378
 379
 380  //  // propagate to the primary vertex
 381 //   if (fInner) {
 382 //     AliExternalTrackParam  *param= new AliExternalTrackParam(*fInner);
 383 //     Double_t pos[]={0.,0.,0.};
 384 //     Double_t sigma[]={.1,.1,.1}; //TODO: is this correct?
 385 //     AliESDVertex vertex(pos,sigma);
 386 //     if (param->PropagateToVertex(&vertex))
 387 //       fPrimary=new AliExternalTrackParam(*param);
 388 //     delete param;
 389 //     // for better comparison one does not want to have alpha changed...
 390 //     if (fPrimary) if (!fPrimary->Rotate(fInner->GetAlpha())) {
 391 //       delete fPrimary;
 392 //       fPrimary=0;
 393 //     }
 394 //   }
 395
 396   delete track;
 397 }
 398
 399
 400
 401
 402 void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
 403                                TrackType type) {
 404   TLinearFitter fy(1);
 405   TLinearFitter fz(1);
 406   fy.StoreData(kFALSE);
 407   fz.StoreData(kFALSE);
 408   switch (type) {
 409   case kLinear:
 410     fy.SetFormula("1 ++ x");
 411     fz.SetFormula("1 ++ x");
 412     break;
 413   case kQuadratic:
 414     fy.SetFormula("1 ++ x ++ x*x");
 415     fz.SetFormula("1 ++ x");
 416     break;
 417   case kKalman:
 418   case kRiemann:
 419     break;
 420   }
 421   Double_t xmax=-1.;
 422   Double_t xmin=1000.;
 423   for (Int_t i=0;i<160;++i) {
 424     AliTPCclusterMI *c=track->GetClusterPointer(i);
 425     if (c && RejectCluster(c)) continue;
 426     if (c&&c->GetDetector()==sector) {
 427       Double_t x=c->GetX();
 428       fy.AddPoint(&x,c->GetY());
 429       fz.AddPoint(&x,c->GetZ());
 430       xmax=TMath::Max(xmax,x);
 431       xmin=TMath::Min(xmin,x);
 432     }
 433   }
 434   fy.Eval();
 435   fz.Eval();
 436   Double_t a[3]={fy.GetParameter(0),
 437                  fy.GetParameter(1),
 438                  type==kQuadratic?fy.GetParameter(2):0.};
 439   Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
 440                   fy.GetCovarianceMatrixElement(1,0),
 441                   fy.GetCovarianceMatrixElement(1,1),
 442                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0.,
 443                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0.,
 444                   type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.};
 445   for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters;
 446   Double_t b[2]={fz.GetParameter(0),
 447                  fz.GetParameter(1)};
 448   Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
 449                   fz.GetCovarianceMatrixElement(1,0),
 450                   fz.GetCovarianceMatrixElement(1,1)};
 451   for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters;
 452   Double_t p[5];
 453   Double_t c[15];
 454   Double_t alpha=track->GetAlpha();
 455   Quadratic2Helix(a,ca,b,cb,0.,p,c);
 456   fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
 457   Quadratic2Helix(a,ca,b,cb,xmin,p,c);
 458   fInner=new AliExternalTrackParam(xmin,alpha,p,c);
 459   Quadratic2Helix(a,ca,b,cb,xmax,p,c);
 460   fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
 461 }
 462
 463 void AliTPCTracklet::Quadratic2Helix(Double_t *a,Double_t *ca,
 464                                      Double_t *b,Double_t *cb,
 465                                      Double_t x0,
 466                                      Double_t *p,Double_t *c) {
 467   // y(x)=a[0]+a[1]*x+a[2]*x^2
 468   // z(x)=b[0]+b[1]*x
 469   // parametrises the corosponding helix at x0
 470
 471   // get the polynoms at x0
 472   Double_t a0=x0*x0*a[2] + x0*a[1] + a[0];
 473   Double_t a1=2.*x0*a[2] +     a[1];
 474   Double_t a2=      a[2];
 475   Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5])));
 476   Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5]));
 477   Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]);
 478   Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]);
 479   Double_t ca21=ca[3]+x0*2.*ca[5];
 480   Double_t ca22=ca[5];
 481
 482   Double_t b0=x0*b[1] + b[0];
 483   Double_t b1=   b[1];
 484   Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]);
 485   Double_t cb10=cb[1]+x0*cb[2];
 486   Double_t cb11=cb[2];
 487
 488   // transform to helix parameters
 489   Double_t f   =1.+a1*a1;
 490   Double_t f2  =f*f;
 491   Double_t fi  =1./f;
 492   Double_t fi12=TMath::Sqrt(fi);
 493   Double_t fi32=fi*fi12;
 494   Double_t fi2 =fi*fi;
 495   Double_t fi52=fi2*fi12;
 496   Double_t fi3 =fi2*fi;
 497   Double_t fi5 =fi2*fi3;
 498
 499   Double_t xyz[3]={0.}; // TODO...
 500   Double_t fc=1./(GetBz(xyz)*kB2C);
 501
 502   p[0]=a0;            // y0
 503   p[1]=b0;            // z0
 504   p[2]=a1*fi12;       // snp
 505   p[3]=b1;            // tgl
 506   p[4]=2.*a2*fi32*fc; // 1/pt
 507
 508   c[0] =ca00;      //  y0-y0
 509   c[1] =0.;        //  z0-y0
 510   c[2] =cb00;      //  z0-z0
 511   c[3] =ca10*fi32; // snp-y0
 512   c[4] =0.;        // snp-z0
 513   c[5] =ca11*fi3;  // snp-snp
 514   c[6] =0.;        // tgl-y0
 515   c[7] =cb10;      // tgl-z0
 516   c[8] =0.;        // tgl-snp
 517   c[9] =cb11;      // tgl-tgl
 518   c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc;  // 1/pt-y0
 519   c[11]=0.;                                 // 1/pt-z0
 520   c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp
 521   c[13]=0.;                                 // 1/pt-tgl
 522   c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5
 523     *fc*fc;        // 1/pt-1/pt
 524 }
 525
 526
 527 void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) {
 528   TLinearFitter fy(2);
 529   fy.StoreData(kFALSE);
 530   fy.SetFormula("hyp2");
 531   Double_t xmax=-1.;
 532   Double_t xmin=1000.;
 533   for (Int_t i=0;i<160;++i) {
 534     AliTPCclusterMI *c=track->GetClusterPointer(i);
 535     if (c && RejectCluster(c)) continue;
 536     if (c&&c->GetDetector()==sector) {
 537       Double_t x=c->GetX();
 538       Double_t y=c->GetY();
 539       Double_t xy[2]={x,y};
 540       Double_t r=x*x+y*y;
 541       Double_t errx=1.,erry=1.;//TODO!
 542       Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry);
 543       err=1.;
 544       fy.AddPoint(xy,r,err);
 545       xmax=TMath::Max(xmax,x);
 546       xmin=TMath::Min(xmin,x);
 547     }
 548   }
 549   fy.Eval();
 550   Double_t a[3]={fy.GetParameter(0),
 551                  fy.GetParameter(1),
 552                  fy.GetParameter(2)};
 553   Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
 554                   fy.GetCovarianceMatrixElement(1,0),
 555                   fy.GetCovarianceMatrixElement(1,1),
 556                   fy.GetCovarianceMatrixElement(2,0),
 557                   fy.GetCovarianceMatrixElement(2,1),
 558                   fy.GetCovarianceMatrixElement(2,2)};
 559
 560   TLinearFitter fz(1);
 561   fz.StoreData(kFALSE);
 562   fz.SetFormula("hyp1");
 563   Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
 564   Double_t oldx=0.;
 565   Double_t oldy=R;
 566   Double_t phi=0.;
 567   for (Int_t i=0;i<160;++i) {
 568     AliTPCclusterMI *c=track->GetClusterPointer(i);
 569     if (c && RejectCluster(c)) continue;
 570     if (c&&c->GetDetector()==sector) {
 571       Double_t x=c->GetX();
 572       Double_t y=c->GetY();
 573       Double_t dx=x-oldx;
 574       Double_t dy=y-oldy;
 575       phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R));
 576       Double_t err=1.;
 577       fz.AddPoint(&phi,c->GetZ(),err);
 578       oldx=x;
 579       oldy=y;
 580     }
 581   }
 582   fz.Eval();
 583   Double_t b[2]={fz.GetParameter(0),
 584                  fz.GetParameter(1)};
 585   Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
 586                   fz.GetCovarianceMatrixElement(1,0),
 587                   fz.GetCovarianceMatrixElement(1,1)};
 588
 589   Double_t p[5];
 590   Double_t c[15];
 591   Double_t alpha=track->GetAlpha();
 592   if (Riemann2Helix(a,ca,b,cb,0.,p,c))
 593     fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
 594   if (Riemann2Helix(a,ca,b,cb,xmin,p,c))
 595     fInner=new AliExternalTrackParam(xmin,alpha,p,c);
 596   if (Riemann2Helix(a,ca,b,cb,xmax,p,c))
 597     fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
 598 }
 599
 600 Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t */*ca*/,
 601                                      Double_t *b,Double_t */*cb*/,
 602                                      Double_t x0,
 603                                      Double_t *p,Double_t *c) {
 604   //TODO: signs!
 605
 606   Double_t xr0=.5*a[1];
 607   Double_t yr0=.5*a[2];
 608   Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
 609   Double_t dx=x0-xr0;
 610   if (dx*dx>=R*R) return kFALSE;
 611   Double_t dy=TMath::Sqrt(R*R-dx*dx); //sign!!
 612   if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
 613     dy=-dy;
 614   Double_t y0=yr0+dy;
 615   Double_t tgp=-dx/dy; //TODO: dy!=0
 616   Double_t z0=b[0]+TMath::ATan(tgp)*b[1];
 617   Double_t xyz[3]={x0,y0,z0};
 618   Double_t fc=1./(GetBz(xyz)*kB2C);
 619   fc=1;
 620   p[0]=y0;  // y0
 621   p[1]=z0; // z0
 622   p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp
 623   p[3]=b[1];       // tgl
 624   p[4]=1./R*fc;    // 1/pt
 625
 626   c[0] =0.;      //  y0-y0
 627   c[1] =0.;        //  z0-y0
 628   c[2] =0.;      //  z0-z0
 629   c[3] =0.; // snp-y0
 630   c[4] =0.;        // snp-z0
 631   c[5] =0.;  // snp-snp
 632   c[6] =0.;        // tgl-y0
 633   c[7] =0.;      // tgl-z0
 634   c[8] =0.;        // tgl-snp
 635   c[9] =0.;      // tgl-tgl
 636   c[10]=0.;  // 1/pt-y0
 637   c[11]=0.;                                 // 1/pt-z0
 638   c[12]=0.; // 1/pt-snp
 639   c[13]=0.;                                 // 1/pt-tgl
 640   c[14]=0.;        // 1/pt-1/pt
 641
 642   return kTRUE;
 643 }
 644
 645 TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
 646                                           TrackType type,
 647                                           Bool_t storeClusters,
 648                                           Int_t minClusters,
 649                                           Int_t maxTracklets) {
 650 // The tracklet factory: It creates several tracklets out of a track. They
 651 // are created for sectors that fullfill the constraint of having enough
 652 // clusters inside. Futhermore you can specify the maximum amount of
 653 // tracklets that are to be created.
 654 // The tracklets appear in a sorted fashion, beginning with those having the
 655 // most clusters.
 656
 657   Int_t sectors[72]={0};
 658   for (Int_t i=0;i<160;++i) {
 659     AliTPCclusterMI *c=track->GetClusterPointer(i);
 660     if (c && RejectCluster(c)) continue;
 661     if (c)
 662       ++sectors[c->GetDetector()];
 663   }
 664   Int_t indices[72];
 665   TMath::Sort(72,sectors,indices);
 666   TObjArray tracklets;
 667   if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
 668   for (Int_t i=0;i<maxTracklets&&sectors[indices[i]]>=minClusters;++i) {
 669     tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
 670   }
 671   return tracklets;
 672 }
 673
 674 TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &/*clusters*/,
 675                                           TrackType /*type*/,
 676                                           Bool_t /*storeClusters*/,
 677                                           Int_t /*minClusters*/,
 678                                           Int_t /*maxTracklets*/) {
 679   // TODO!
 680
 681   TObjArray tracklets;
 682   return tracklets;
 683 }
 684
 685 Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
 686                                         const AliTPCTracklet &t2,
 687                                         AliExternalTrackParam *&t1m,
 688                                         AliExternalTrackParam *&t2m) {
 689   // This function propagates two Tracklets to a common x-coordinate. This
 690   // x is dermined as the one that is in the middle of the two tracklets (they
 691   // are assumed to live on two distinct x-intervalls).
 692   // The inner parametrisation of the outer Tracklet and the outer
 693   // parametrisation of the inner Tracklet are used and propagated to this
 694   // common x. This result is saved not inside the Tracklets but two new
 695   // ExternalTrackParams are created (that means you might want to delete
 696   // them).
 697   // In the case that the alpha angles of the Tracklets differ both angles
 698   // are tried out for this propagation.
 699   // In case of any failure kFALSE is returned, no AliExternalTrackParam
 700   // is created und the pointers are set to 0.
 701
 702   if (t1.GetInner() && t1.GetOuter() &&
 703       t2.GetInner() && t2.GetOuter()) {
 704     if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) {
 705       t1m=new AliExternalTrackParam(*t1.GetOuter());
 706       t2m=new AliExternalTrackParam(*t2.GetInner());
 707     }
 708     else {
 709       t1m=new AliExternalTrackParam(*t1.GetInner());
 710       t2m=new AliExternalTrackParam(*t2.GetOuter());
 711     }
 712     Double_t mx=.5*(t1m->GetX()+t2m->GetX());
 713     Double_t b1,b2;
 714     Double_t xyz[3];
 715     t1m->GetXYZ(xyz);
 716     b1=GetBz(xyz);
 717     t2m->GetXYZ(xyz);
 718     b2=GetBz(xyz);
 719     if (t1m->Rotate(t2m->GetAlpha())
 720         && t1m->PropagateTo(mx,b1)
 721         && t2m->PropagateTo(mx,b2));
 722     else
 723       if (t2m->Rotate(t1m->GetAlpha())
 724           && t1m->PropagateTo(mx,b1)
 725           && t2m->PropagateTo(mx,b2));
 726       else {
 727         delete t1m;
 728         delete t2m;
 729         t1m=t2m=0;
 730       }
 731   }
 732   else {
 733     t1m=t2m=0;
 734   }
 735   return t1m&&t2m;
 736 }
 737
 738 double AliTPCTracklet::GetBz(Double_t *xyz) {
 739   if (AliTracker::UniformField())
 740     return AliTracker::GetBz();
 741   else
 742     return AliTracker::GetBz(xyz);
 743 }
 744
 745 void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c,
 746                               Double_t *x) {
 747   // This function generates a n-dimensional random variable x with mean
 748   // p and covariance c.
 749   // That is done using the cholesky decomposition of the covariance matrix,
 750   // Begin_Latex C=U^{t} U End_Latex, with Begin_Latex U End_Latex being an
 751   // upper triangular matrix. Given a vector v of iid gausian random variables
 752   // with variance 1 one obtains the asked result as: Begin_Latex x=U^t v
 753   // End_Latex.
 754   // c is expected to be in a lower triangular format:
 755   // c[0]
 756   // c[1] c[2]
 757   // c[3] c[4] c[5]
 758   // etc.
 759   static TRandom3 random;
 760   Double_t *c2= new Double_t[ndim*ndim];
 761   Int_t k=0;
 762   for (Int_t i=0;i<ndim;++i)
 763     for (Int_t j=0;j<=i;++j)
 764       c2[i*ndim+j]=c2[j*ndim+i]=c[k++];
 765   TMatrixDSym cm(ndim,c2);
 766   delete[] c2;
 767   TDecompChol chol(cm);
 768   chol.Decompose();
 769   const TVectorD pv(ndim);
 770   const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p));
 771   TVectorD xv(ndim);
 772   xv.Use(ndim,x);
 773   for (Int_t i=0;i<ndim;++i)
 774     xv[i]=random.Gaus();
 775   TMatrixD L=chol.GetU();
 776   L.T();
 777   xv=L*xv+pv;
 778 }
 779
 780 TEllipse AliTPCTracklet::ErrorEllipse(Double_t x,Double_t y,
 781                                       Double_t sx,Double_t sy,Double_t sxy) {
 782   /* Begin_Latex
 783      r_{1,2}=1/2 (a+c#pm#sqrt{(a-c)^{2}+(2b)^{2}})
 784   End_Latex */
 785   Double_t det1=1./(sx*sy-sxy*sxy);
 786   Double_t a=sy*det1;
 787   Double_t b=-sxy*det1;
 788   Double_t c=sx*det1;
 789   Double_t d=c-a;
 790   Double_t s=TMath::Sqrt(d*d+4.*b*b);
 791   Double_t r1=TMath::Sqrt(.5*(a+c-s));
 792   Double_t r2=TMath::Sqrt(.5*(a+c+s));
 793   Double_t alpha=.5*TMath::ATan2(2.*b,d);
 794   return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg());
 795 }
 796
 797 void AliTPCTracklet::Test(const char* filename) {
 798   /*
 799     aliroot
 800     AliTPCTracklet::Test("");
 801     TFile f("AliTPCTrackletDebug.root");
 802     TTree *t=f.Get("AliTPCTrackletDebug");
 803     t->Draw("p0:p4");
 804     TEllipse e=AliTPCTracklet::ErrorEllipse(0.,0.,4.,1.,1.8);
 805     e.Draw();
 806  */
 807   TTreeSRedirector ds(filename);
 808   Double_t p[5]={0.};
 809   Double_t c[15]={4.,
 810                   0.,4.,
 811                   0.,0.,9.,
 812                   0.,0.,0.,16.,
 813                   1.8,0.,0.,0.,1.};
 814   for (Int_t i=0;i<10000;++i) {
 815     Double_t x[5];
 816     RandomND(5,p,c,x);
 817     ds<<"AliTPCTrackletDebug"
 818       <<"p0="<<x[0]
 819       <<"p1="<<x[1]
 820       <<"p2="<<x[2]
 821       <<"p3="<<x[3]
 822       <<"p4="<<x[4]
 823       <<"\n";
 824   }
 825
 826   /*
 827   Double_t b;
 828   Double_t x=0.;
 829   Double_t alpha=0.;
 830   Double_t param[5]={0.};
 831   Double_t covar[15]={1.,
 832                       0.,1.,
 833                       0.,0.,1.,
 834                       0.,0.,0.,1.,
 835                       0.,0.,0.,0.,1.};
 836   AliExternalTrackParam track(x,alpha,param,covar);
 837
 838
 839
 840   for (Int_t i=0;i<points.GetNPoints();++i) {
 841     Double_t x=0.;
 842     Double_t alpha=0.;
 843     Double_t param[5]={0.};
 844     Double_t covar[15]={1.,
 845                         0.,1.,
 846                         0.,0.,1.,
 847                         0.,0.,0.,1.,
 848                         0.,0.,0.,0.,1.};
 849     AliExternalTrackParam track(x,alpha,param,covar);
 850     for (x=90.;x<250.;x+=1.) {
 851       track.PropagateTo(x,b);
 852       AliTPCclusterMI c();
 853     }
 854   }
 855   */
 856 }
 857
 858
 859 Bool_t AliTPCTracklet::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
 860   //
 861   // check the acceptance of cluster
 862   // Cut on edge effects
 863   //
 864   Bool_t isReject = kFALSE;
 865   Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
 866   Float_t dist  = edgeY - TMath::Abs(cl->GetY());
 867   if (param)  dist  = edgeY - TMath::Abs(param->GetY());
 868   if (dist<fgEdgeCutY) isReject=kTRUE;
 869   if (cl->GetType()<0) isReject=kTRUE;
 870   return isReject;
 871 }