TPC/AliTPCtracker.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 /* $Id$ */
  17
  18 //-------------------------------------------------------
  19 //          Implementation of the TPC tracker
  20 //
  21 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  22 //-------------------------------------------------------
  23 #include <TObjArray.h>
  24 #include <TError.h>
  25 #include <TFile.h>
  26 #include <TTree.h>
  27
  28 #include "AliESD.h"
  29
  30 #include "AliTPCtracker.h"
  31 #include "AliTPCcluster.h"
  32 #include "AliTPCParam.h"
  33 #include "AliClusters.h"
  34
  35 ClassImp(AliTPCtracker)
  36
  37 //_____________________________________________________________________________
  38 AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
  39 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
  40 {
  41   //---------------------------------------------------------------------
  42   // The main TPC tracker constructor
  43   //---------------------------------------------------------------------
  44   fInnerSec=new AliTPCSector[fkNIS];
  45   fOuterSec=new AliTPCSector[fkNOS];
  46
  47   Int_t i;
  48   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
  49   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
  50
  51   fParam = (AliTPCParam*) par;
  52   fSeeds=0;
  53
  54   // [SR 17.03.2003]
  55
  56   fBarrelFile = 0;
  57   fBarrelTree = 0;
  58   fBarrelArray = 0;
  59   fBarrelTrack = 0;
  60
  61 }
  62
  63 //_____________________________________________________________________________
  64 AliTPCtracker::~AliTPCtracker() {
  65   //------------------------------------------------------------------
  66   // TPC tracker destructor
  67   //------------------------------------------------------------------
  68   delete[] fInnerSec;
  69   delete[] fOuterSec;
  70   if (fSeeds) {
  71     fSeeds->Delete();
  72     delete fSeeds;
  73   }
  74
  75   // [SR, 01.04.2003]
  76   if (fBarrelFile) {
  77     fBarrelFile->Close();
  78     delete fBarrelFile;
  79   }
  80 }
  81
  82 //_____________________________________________________________________________
  83
  84 void AliTPCtracker::SetBarrelTree(const char *mode) {
  85   //
  86   // Creates a tree for BarrelTracks
  87   // mode = "back" or "refit"
  88   //
  89   // [SR, 01.04.2003]
  90   //
  91
  92   if (!IsStoringBarrel()) return;
  93
  94   TDirectory *sav = gDirectory;
  95   if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
  96
  97   char buff[40];
  98   sprintf(buff, "BarrelTPC_%d_%s", GetEventNumber(), mode);
  99
 100   fBarrelFile->cd();
 101   fBarrelTree = new TTree(buff, "Barrel TPC tracks");
 102
 103   if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", 4);
 104   for(Int_t i=0; i<4; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
 105
 106   fBarrelTree->Branch("tracks", &fBarrelArray);
 107
 108   sav->cd();
 109 }
 110 //_____________________________________________________________________________
 111
 112 void AliTPCtracker::StoreBarrelTrack(AliTPCtrack *ps, Int_t refPlane, Int_t isIn) {
 113   //
 114   // Stores Track at a given reference plane
 115   //
 116   // refPlane: 1-4
 117   // isIn: 1 - backward, 2 - refit
 118   //
 119
 120   if (!IsStoringBarrel()) return;
 121   if (refPlane < 0 || refPlane > 4) return;
 122   if (isIn > 2) return;
 123
 124   static Int_t nClusters;
 125   static Int_t nWrong;
 126   static Double_t chi2;
 127   static Int_t index;
 128
 129   Int_t newClusters, newWrong;
 130   Double_t newChi2;
 131
 132   if ( (refPlane == 1 && isIn == kTrackBack) ||
 133        (refPlane == 4 && isIn == kTrackRefit) ) {
 134
 135     fBarrelArray->Clear();
 136     nClusters = nWrong = 0;
 137     chi2 = 0.0;
 138     index = 0;
 139   }
 140
 141   // propagate
 142   Double_t refX = 0;
 143   if (refPlane == 1) refX = fParam->GetInnerRadiusLow();
 144   if (refPlane == 2) refX = fParam->GetInnerRadiusUp();
 145   if (refPlane == 3) refX = fParam->GetOuterRadiusLow();
 146   if (refPlane == 4) refX = fParam->GetOuterRadiusUp();
 147
 148   ps->PropagateTo(refX);
 149
 150   fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
 151   ps->GetBarrelTrack(fBarrelTrack);
 152
 153   newClusters = ps->GetNumberOfClusters() - nClusters;
 154   newWrong = ps->GetNWrong() - nWrong;
 155   newChi2 = ps->GetChi2() - chi2;
 156
 157   nClusters =  ps->GetNumberOfClusters();
 158   nWrong = ps->GetNWrong();
 159   chi2 = ps->GetChi2();
 160
 161   fBarrelTrack->SetNClusters(newClusters, newChi2);
 162   fBarrelTrack->SetNWrongClusters(newWrong);
 163   fBarrelTrack->SetRefPlane(refPlane, isIn);
 164 }
 165
 166
 167 //_____________________________________________________________________________
 168 Double_t f1(Double_t x1,Double_t y1,
 169                    Double_t x2,Double_t y2,
 170                    Double_t x3,Double_t y3)
 171 {
 172   //-----------------------------------------------------------------
 173   // Initial approximation of the track curvature
 174   //-----------------------------------------------------------------
 175   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 176   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 177                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 178   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 179                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 180
 181   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 182
 183   return -xr*yr/sqrt(xr*xr+yr*yr);
 184 }
 185
 186
 187 //_____________________________________________________________________________
 188 Double_t f2(Double_t x1,Double_t y1,
 189                    Double_t x2,Double_t y2,
 190                    Double_t x3,Double_t y3)
 191 {
 192   //-----------------------------------------------------------------
 193   // Initial approximation of the track curvature times center of curvature
 194   //-----------------------------------------------------------------
 195   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 196   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 197                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 198   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 199                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 200
 201   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 202
 203   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 204 }
 205
 206 //_____________________________________________________________________________
 207 Double_t f3(Double_t x1,Double_t y1,
 208                    Double_t x2,Double_t y2,
 209                    Double_t z1,Double_t z2)
 210 {
 211   //-----------------------------------------------------------------
 212   // Initial approximation of the tangent of the track dip angle
 213   //-----------------------------------------------------------------
 214   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 215 }
 216
 217 //_____________________________________________________________________________
 218 Int_t AliTPCtracker::LoadClusters(TTree *cTree) {
 219   //-----------------------------------------------------------------
 220   // This function loads TPC clusters.
 221   //-----------------------------------------------------------------
 222   TBranch *branch=cTree->GetBranch("Segment");
 223   if (!branch) {
 224      Error("LoadClusters","Can't get the branch !");
 225      return 1;
 226   }
 227
 228   AliClusters carray, *addr=&carray;
 229   carray.SetClass("AliTPCcluster");
 230   carray.SetArray(0);
 231   branch->SetAddress(&addr);
 232
 233   Int_t nentr=(Int_t)cTree->GetEntries();
 234
 235   for (Int_t i=0; i<nentr; i++) {
 236       cTree->GetEvent(i);
 237
 238       Int_t ncl=carray.GetArray()->GetEntriesFast();
 239
 240       Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
 241       Int_t id=carray.GetID();
 242       if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
 243          Fatal("LoadClusters","Wrong index !");
 244       }
 245       Int_t outindex = 2*fkNIS*nir;
 246       if (id<outindex) {
 247          Int_t sec = id/nir;
 248          Int_t row = id - sec*nir;
 249          sec %= fkNIS;
 250          AliTPCRow &padrow=fInnerSec[sec][row];
 251          while (ncl--) {
 252            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
 253            padrow.InsertCluster(c,sec,row);
 254          }
 255       } else {
 256          id -= outindex;
 257          Int_t sec = id/nor;
 258          Int_t row = id - sec*nor;
 259          sec %= fkNOS;
 260          AliTPCRow &padrow=fOuterSec[sec][row];
 261          while (ncl--) {
 262            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
 263            padrow.InsertCluster(c,sec+fkNIS,row);
 264          }
 265       }
 266       carray.GetArray()->Clear();
 267   }
 268
 269   return 0;
 270 }
 271
 272 //_____________________________________________________________________________
 273 void AliTPCtracker::UnloadClusters() {
 274   //-----------------------------------------------------------------
 275   // This function unloads TPC clusters.
 276   //-----------------------------------------------------------------
 277   Int_t i;
 278   for (i=0; i<fkNIS; i++) {
 279     Int_t nr=fInnerSec->GetNRows();
 280     for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
 281   }
 282   for (i=0; i<fkNOS; i++) {
 283     Int_t nr=fOuterSec->GetNRows();
 284     for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
 285   }
 286 }
 287
 288 //_____________________________________________________________________________
 289 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
 290   //-----------------------------------------------------------------
 291   // This function tries to find a track prolongation.
 292   //-----------------------------------------------------------------
 293   Double_t xt=t.GetX();
 294   const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
 295                     Int_t(0.5*fSectors->GetNRows());
 296   Int_t tryAgain=kSKIP;
 297
 298   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
 299   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 300   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 301   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 302
 303   Int_t nrows=fSectors->GetRowNumber(xt)-1;
 304   for (Int_t nr=nrows; nr>=rf; nr--) {
 305     Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
 306     if (!t.PropagateTo(x)) return 0;
 307
 308     AliTPCcluster *cl=0;
 309     UInt_t index=0;
 310     Double_t maxchi2=kMaxCHI2;
 311     const AliTPCRow &krow=fSectors[s][nr];
 312     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
 313     Double_t sy2=AliTPCcluster::SigmaY2(t.GetX(),t.GetTgl(),pt);
 314     Double_t sz2=AliTPCcluster::SigmaZ2(t.GetX(),t.GetTgl());
 315     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 316
 317     if (road>kMaxROAD) {
 318       if (t.GetNumberOfClusters()>4)
 319          Warning("FindProlongation","Too broad road !");
 320       return 0;
 321     }
 322
 323     if (krow) {
 324       for (Int_t i=krow.Find(y-road); i<krow; i++) {
 325         AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 326         if (c->GetY() > y+road) break;
 327         if (c->IsUsed()) continue;
 328         if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
 329         Double_t chi2=t.GetPredictedChi2(c);
 330         if (chi2 > maxchi2) continue;
 331         maxchi2=chi2;
 332         cl=c;
 333         index=krow.GetIndex(i);
 334       }
 335     }
 336     if (cl) {
 337       Float_t l=fSectors->GetPadPitchWidth();
 338       Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
 339       t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
 340       if (!t.Update(cl,maxchi2,index)) {
 341          if (!tryAgain--) return 0;
 342       } else tryAgain=kSKIP;
 343     } else {
 344       if (tryAgain==0) break;
 345       if (y > ymax) {
 346          s = (s+1) % fN;
 347          if (!t.Rotate(fSectors->GetAlpha())) return 0;
 348       } else if (y <-ymax) {
 349          s = (s-1+fN) % fN;
 350          if (!t.Rotate(-fSectors->GetAlpha())) return 0;
 351       }
 352       tryAgain--;
 353     }
 354   }
 355
 356   return 1;
 357 }
 358 //_____________________________________________________________________________
 359
 360 Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
 361   //
 362   // This function propagates seed inward TPC using old clusters
 363   // from the track.
 364   //
 365   // Sylwester Radomski, GSI
 366   // 26.02.2003
 367   //
 368
 369   // loop over rows
 370
 371   Int_t nRows = fSectors->GetNRows();
 372   for (Int_t iRow = nRows-1; iRow >= 0; iRow--) {
 373
 374     Double_t x = fSectors->GetX(iRow);
 375     if (!seed->PropagateTo(x)) return 0;
 376
 377     // try to find an assigned cluster in this row
 378
 379     AliTPCcluster* cluster = NULL;
 380     Int_t idx = -1;
 381     Int_t sec = -1;
 382     for (Int_t iCluster = 0; iCluster < track->GetNumberOfClusters(); iCluster++) {
 383       idx = track->GetClusterIndex(iCluster);
 384       sec = (idx&0xff000000)>>24;
 385       Int_t row = (idx&0x00ff0000)>>16;
 386       if (((fSectors == fInnerSec) && (sec >= fkNIS)) ||
 387           ((fSectors == fOuterSec) && (sec < fkNIS))) continue;
 388       if (row == iRow) {
 389         cluster = (AliTPCcluster*) GetCluster(idx);
 390         break;
 391       }
 392     }
 393
 394     // update the track seed with the found cluster
 395
 396     if (cluster) {
 397       Double_t dAlpha = fParam->GetAngle(sec) - seed->GetAlpha();
 398       if (TMath::Abs(dAlpha) > 0.0001) {
 399         if (!seed->Rotate(dAlpha)) return 0;
 400         if (!seed->PropagateTo(x)) return 0;
 401       }
 402
 403       seed->Update(cluster, seed->GetPredictedChi2(cluster), idx);
 404     }
 405   }
 406
 407   return 1;
 408 }
 409
 410 //_____________________________________________________________________________
 411 Int_t AliTPCtracker::FollowBackProlongation
 412 (AliTPCseed& seed, const AliTPCtrack &track) {
 413   //-----------------------------------------------------------------
 414   // This function propagates tracks back through the TPC
 415   //-----------------------------------------------------------------
 416   Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
 417   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 418   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 419   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 420
 421   Int_t idx=-1, sec=-1, row=-1;
 422   Int_t nc=seed.GetNumber();
 423
 424   if (nc--) {
 425      idx=track.GetClusterIndex(nc);
 426      sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 427   }
 428   if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 429   else { if (sec <  fkNIS) row=-1; }
 430
 431   Int_t nr=fSectors->GetNRows();
 432   for (Int_t i=0; i<nr; i++) {
 433     Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
 434     Double_t y=seed.GetYat(x);
 435
 436     if (y > ymax) {
 437        s = (s+1) % fN;
 438        if (!seed.Rotate(fSectors->GetAlpha())) return 0;
 439     } else if (y <-ymax) {
 440        s = (s-1+fN) % fN;
 441        if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
 442     }
 443
 444     if (!seed.PropagateTo(x)) return 0;
 445
 446     AliTPCcluster *cl=0;
 447     Int_t index=0;
 448     Double_t maxchi2=kMaxCHI2;
 449     Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
 450     Double_t sy2=AliTPCcluster::SigmaY2(seed.GetX(),seed.GetTgl(),pt);
 451     Double_t sz2=AliTPCcluster::SigmaZ2(seed.GetX(),seed.GetTgl());
 452     Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
 453     if (road>kMaxROAD) {
 454       Warning("FollowBackProlongation","Too broad road !");
 455       return 0;
 456     }
 457
 458     Int_t accepted=seed.GetNumberOfClusters();
 459     if (row==i) {
 460        //try to accept already found cluster
 461        AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
 462        Double_t chi2;
 463        if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
 464           index=idx; cl=c; maxchi2=chi2;
 465        } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
 466
 467        if (nc--) {
 468           idx=track.GetClusterIndex(nc);
 469           sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 470        }
 471        if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 472        else { if (sec < fkNIS) row=-1; }
 473
 474     }
 475     if (!cl) {
 476        //try to fill the gap
 477        const AliTPCRow &krow=fSectors[s][i];
 478        if (accepted>27)
 479        if (krow) {
 480           for (Int_t i=krow.Find(y-road); i<krow; i++) {
 481             AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 482             if (c->GetY() > y+road) break;
 483             if (c->IsUsed()) continue;
 484          if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
 485             Double_t chi2=seed.GetPredictedChi2(c);
 486             if (chi2 > maxchi2) continue;
 487             maxchi2=chi2;
 488             cl=c;
 489             index=krow.GetIndex(i);
 490           }
 491        }
 492     }
 493
 494     if (cl) {
 495       Float_t l=fSectors->GetPadPitchWidth();
 496       Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
 497       seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
 498       seed.Update(cl,maxchi2,index);
 499     }
 500
 501   }
 502
 503   seed.SetNumber(nc);
 504
 505   return 1;
 506 }
 507
 508 //_____________________________________________________________________________
 509 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
 510   //-----------------------------------------------------------------
 511   // This function creates track seeds.
 512   //-----------------------------------------------------------------
 513   Double_t x[5], c[15];
 514
 515   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
 516   Double_t cs=cos(alpha), sn=sin(alpha);
 517
 518   Double_t x1 =fSectors->GetX(i1);
 519   Double_t xx2=fSectors->GetX(i2);
 520
 521   for (Int_t ns=0; ns<fN; ns++) {
 522     Int_t nl=fSectors[(ns-1+fN)%fN][i2];
 523     Int_t nm=fSectors[ns][i2];
 524     Int_t nu=fSectors[(ns+1)%fN][i2];
 525     const AliTPCRow& kr1=fSectors[ns][i1];
 526     for (Int_t is=0; is < kr1; is++) {
 527       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
 528       for (Int_t js=0; js < nl+nm+nu; js++) {
 529         const AliTPCcluster *kcl;
 530         Double_t x2,   y2,   z2;
 531         Double_t x3=GetX(), y3=GetY(), z3=GetZ();
 532
 533         if (js<nl) {
 534           const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
 535           kcl=kr2[js];
 536           y2=kcl->GetY(); z2=kcl->GetZ();
 537           x2= xx2*cs+y2*sn;
 538           y2=-xx2*sn+y2*cs;
 539         } else
 540           if (js<nl+nm) {
 541             const AliTPCRow& kr2=fSectors[ns][i2];
 542             kcl=kr2[js-nl];
 543             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
 544           } else {
 545             const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
 546             kcl=kr2[js-nl-nm];
 547             y2=kcl->GetY(); z2=kcl->GetZ();
 548             x2=xx2*cs-y2*sn;
 549             y2=xx2*sn+y2*cs;
 550           }
 551
 552         Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
 553         if (TMath::Abs(zz-z2)>5.) continue;
 554
 555         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
 556         if (d==0.) {
 557            Warning("MakeSeeds","Straight seed !");
 558            continue;
 559         }
 560         x[0]=y1;
 561         x[1]=z1;
 562         x[4]=f1(x1,y1,x2,y2,x3,y3);
 563         if (TMath::Abs(x[4]) >= 0.0066) continue;
 564         x[2]=f2(x1,y1,x2,y2,x3,y3);
 565         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
 566         x[3]=f3(x1,y1,x2,y2,z1,z2);
 567         if (TMath::Abs(x[3]) > 1.2) continue;
 568         Double_t a=asin(x[2]);
 569         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
 570         if (TMath::Abs(zv-z3)>10.) continue;
 571
 572         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
 573         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
 574         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
 575         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
 576
 577         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
 578         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
 579         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
 580         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
 581         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
 582         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
 583         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
 584         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
 585         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
 586         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
 587
 588         c[0]=sy1;
 589         c[1]=0.;       c[2]=sz1;
 590         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
 591         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
 592                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
 593         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
 594         c[13]=f30*sy1*f40+f32*sy2*f42;
 595         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
 596
 597         UInt_t index=kr1.GetIndex(is);
 598         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
 599         Float_t l=fSectors->GetPadPitchWidth();
 600         track->SetSampledEdx(kr1[is]->GetQ()/l,0);
 601
 602         Int_t rc=FollowProlongation(*track, i2);
 603         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
 604         else fSeeds->AddLast(track);
 605       }
 606     }
 607   }
 608 }
 609
 610 //_____________________________________________________________________________
 611 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
 612   //-----------------------------------------------------------------
 613   // This function reades track seeds.
 614   //-----------------------------------------------------------------
 615   TDirectory *savedir=gDirectory;
 616
 617   TFile *in=(TFile*)inp;
 618   if (!in->IsOpen()) {
 619      Error("ReadSeeds","Input file has not been open !");
 620      return 1;
 621   }
 622
 623   in->cd();
 624   TTree *seedTree=(TTree*)in->Get("Seeds");
 625   if (!seedTree) {
 626      Error("ReadSeeds","Can't get a tree with track seeds !");
 627      return 2;
 628   }
 629   AliTPCtrack *seed=new AliTPCtrack;
 630   seedTree->SetBranchAddress("tracks",&seed);
 631
 632   if (fSeeds==0) fSeeds=new TObjArray(15000);
 633
 634   Int_t n=(Int_t)seedTree->GetEntries();
 635   for (Int_t i=0; i<n; i++) {
 636      seedTree->GetEvent(i);
 637      fSeeds->AddLast(new AliTPCseed(*seed));
 638   }
 639
 640   delete seed;
 641
 642   delete seedTree; //Thanks to Mariana Bondila
 643
 644   savedir->cd();
 645
 646   return 0;
 647 }
 648
 649 //_____________________________________________________________________________
 650 Int_t AliTPCtracker::Clusters2Tracks(AliESD *event) {
 651   //-----------------------------------------------------------------
 652   // This is a track finder.
 653   // The clusters must be already loaded !
 654   //-----------------------------------------------------------------
 655
 656   //find track seeds
 657   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 658   Int_t nrows=nlow+nup;
 659   if (fSeeds==0) {
 660      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
 661      fSectors=fOuterSec; fN=fkNOS;
 662      fSeeds=new TObjArray(15000);
 663      MakeSeeds(nup-1, nup-1-gap);
 664      MakeSeeds(nup-1-shift, nup-1-shift-gap);
 665   }
 666   fSeeds->Sort();
 667
 668   Int_t nseed=fSeeds->GetEntriesFast();
 669   for (Int_t i=0; i<nseed; i++) {
 670     //tracking in the outer sectors
 671     fSectors=fOuterSec; fN=fkNOS;
 672
 673     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 674     if (!FollowProlongation(t)) {
 675        delete fSeeds->RemoveAt(i);
 676        continue;
 677     }
 678
 679     //tracking in the inner sectors
 680     fSectors=fInnerSec; fN=fkNIS;
 681
 682     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
 683     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 684     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 685     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
 686
 687     alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
 688
 689     if (t.Rotate(alpha)) {
 690       if (FollowProlongation(t)) {
 691         if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 692           t.CookdEdx();
 693           CookLabel(pt,0.1); //For comparison only
 694           pt->PropagateTo(fParam->GetInnerRadiusLow());
 695           AliESDtrack iotrack;
 696           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 697
 698           event->AddTrack(&iotrack);
 699
 700           UseClusters(&t);
 701         }
 702       }
 703     }
 704     delete fSeeds->RemoveAt(i);
 705   }
 706
 707   Info("Clusters2Tracks","Number of found tracks : %d",
 708        event->GetNumberOfTracks());
 709
 710   fSeeds->Clear(); delete fSeeds; fSeeds=0;
 711
 712   return 0;
 713 }
 714
 715 //_____________________________________________________________________________
 716 Int_t AliTPCtracker::RefitInward(AliESD* event) {
 717   //
 718   // The function propagates tracks throught TPC inward
 719   // using already associated clusters.
 720   // The clusters must be already loaded !
 721   //
 722
 723   Int_t nTracks = event->GetNumberOfTracks();
 724   Int_t nRefited = 0;
 725
 726   for (Int_t i = 0; i < nTracks; i++) {
 727     AliESDtrack* track = event->GetTrack(i);
 728     ULong_t status = track->GetStatus();
 729
 730     if ( (status & AliESDtrack::kTPCrefit) != 0 ) continue;
 731     if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
 732
 733     if ( (status & AliESDtrack::kTRDout ) != 0 )
 734       if ( (status & AliESDtrack::kTRDrefit ) == 0 ) continue;
 735
 736     AliTPCtrack* tpcTrack = new AliTPCtrack(*track);
 737     AliTPCseed* seed = new AliTPCseed(*tpcTrack);
 738
 739     if ( (status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
 740
 741     fSectors = fOuterSec;
 742
 743     Int_t res = FollowRefitInward(seed, tpcTrack);
 744     UseClusters(seed);
 745     Int_t nc = seed->GetNumberOfClusters();
 746
 747     fSectors = fInnerSec;
 748
 749     res = FollowRefitInward(seed, tpcTrack);
 750     UseClusters(seed, nc);
 751
 752     if (res) {
 753       seed->PropagateTo(fParam->GetInnerRadiusLow());
 754       seed->SetLabel(tpcTrack->GetLabel());
 755       seed->SetdEdx(tpcTrack->GetdEdx());
 756       track->UpdateTrackParams(seed, AliESDtrack::kTPCrefit);
 757       nRefited++;
 758     }
 759
 760     delete seed;
 761     delete tpcTrack;
 762   }
 763
 764   Info("RefitInward","Number of refitted tracks : %d",nRefited);
 765
 766   return 0;
 767 }
 768
 769 Int_t AliTPCtracker::PropagateBack(AliESD *event) {
 770   //-----------------------------------------------------------------
 771   // This function propagates tracks back through the TPC.
 772   // The clusters must be already loaded !
 773   //-----------------------------------------------------------------
 774   Int_t nentr=event->GetNumberOfTracks();
 775   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 776
 777   Int_t ntrk=0;
 778   for (Int_t i=0; i<nentr; i++) {
 779     AliESDtrack *esd=event->GetTrack(i);
 780     ULong_t status=esd->GetStatus();
 781
 782     if ( (status & AliESDtrack::kTPCin ) == 0 ) continue;
 783     if ( (status & AliESDtrack::kTPCout) != 0 ) continue;
 784     if ( (status & AliESDtrack::kITSin) != 0 )
 785        if ( (status & AliESDtrack::kITSout) == 0 ) continue;
 786
 787     const AliTPCtrack t(*esd);
 788     AliTPCseed s(t);
 789
 790     if ( (status & AliESDtrack::kITSout) == 0 ) s.ResetCovariance();
 791
 792     s.ResetNWrong();
 793     s.ResetNRotation();
 794
 795     Int_t nc=t.GetNumberOfClusters();
 796     s.SetNumber(nc); //set number of the cluster to start with
 797
 798     //inner sectors
 799     fSectors=fInnerSec; fN=fkNIS;
 800
 801     Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
 802     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 803     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 804     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
 805     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
 806     alpha-=s.GetAlpha();
 807
 808     if (!s.Rotate(alpha)) continue;
 809     if (!FollowBackProlongation(s,t)) continue;
 810
 811     UseClusters(&s);
 812
 813     //outer sectors
 814     fSectors=fOuterSec; fN=fkNOS;
 815
 816     nc=s.GetNumberOfClusters();
 817
 818     alpha=s.GetAlpha() - fSectors->GetAlphaShift();
 819     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 820     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 821     ns=Int_t(alpha/fSectors->GetAlpha())%fN;
 822
 823     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
 824     alpha-=s.GetAlpha();
 825
 826     if (!s.Rotate(alpha)) continue;
 827     if (!FollowBackProlongation(s,t)) continue;
 828     {
 829     Int_t nrows=fOuterSec->GetNRows()+fInnerSec->GetNRows();
 830     if (s.GetNumberOfClusters() < Int_t(0.4*nrows)) continue;
 831     }
 832     s.PropagateTo(fParam->GetOuterRadiusUp());
 833     s.CookdEdx();
 834     CookLabel(&s,0.1); //For comparison only
 835     UseClusters(&s,nc);
 836     esd->UpdateTrackParams(&s,AliESDtrack::kTPCout);
 837     ntrk++;
 838   }
 839   Info("PropagateBack","Number of back propagated tracks: %d",ntrk);
 840
 841   return 0;
 842 }
 843
 844 //_________________________________________________________________________
 845 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
 846   //--------------------------------------------------------------------
 847   //       Return pointer to a given cluster
 848   //--------------------------------------------------------------------
 849   Int_t sec=(index&0xff000000)>>24;
 850   Int_t row=(index&0x00ff0000)>>16;
 851   Int_t ncl=(index&0x0000ffff)>>00;
 852
 853   const AliTPCcluster *cl=0;
 854
 855   if (sec<fkNIS) {
 856     cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
 857   } else {
 858     sec-=fkNIS;
 859     cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
 860   }
 861
 862   return (AliCluster*)cl;
 863 }
 864
 865 //__________________________________________________________________________
 866 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
 867   //--------------------------------------------------------------------
 868   //This function "cooks" a track label. If label<0, this track is fake.
 869   //--------------------------------------------------------------------
 870   Int_t noc=t->GetNumberOfClusters();
 871   Int_t *lb=new Int_t[noc];
 872   Int_t *mx=new Int_t[noc];
 873   AliCluster **clusters=new AliCluster*[noc];
 874
 875   Int_t i;
 876   for (i=0; i<noc; i++) {
 877      lb[i]=mx[i]=0;
 878      Int_t index=t->GetClusterIndex(i);
 879      clusters[i]=GetCluster(index);
 880   }
 881
 882   Int_t lab=123456789;
 883   for (i=0; i<noc; i++) {
 884     AliCluster *c=clusters[i];
 885     lab=TMath::Abs(c->GetLabel(0));
 886     Int_t j;
 887     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
 888     lb[j]=lab;
 889     (mx[j])++;
 890   }
 891
 892   Int_t max=0;
 893   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
 894
 895   for (i=0; i<noc; i++) {
 896     AliCluster *c=clusters[i];
 897     if (TMath::Abs(c->GetLabel(1)) == lab ||
 898         TMath::Abs(c->GetLabel(2)) == lab ) max++;
 899   }
 900
 901   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
 902
 903   else {
 904      Int_t tail=Int_t(0.10*noc);
 905      max=0;
 906      for (i=1; i<=tail; i++) {
 907        AliCluster *c=clusters[noc-i];
 908        if (lab == TMath::Abs(c->GetLabel(0)) ||
 909            lab == TMath::Abs(c->GetLabel(1)) ||
 910            lab == TMath::Abs(c->GetLabel(2))) max++;
 911      }
 912      if (max < Int_t(0.5*tail)) lab=-lab;
 913   }
 914
 915   t->SetLabel(lab);
 916
 917   delete[] lb;
 918   delete[] mx;
 919   delete[] clusters;
 920 }
 921
 922 //_________________________________________________________________________
 923 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
 924   //-----------------------------------------------------------------------
 925   // Setup inner sector
 926   //-----------------------------------------------------------------------
 927   if (f==0) {
 928      fAlpha=par->GetInnerAngle();
 929      fAlphaShift=par->GetInnerAngleShift();
 930      fPadPitchWidth=par->GetInnerPadPitchWidth();
 931      f1PadPitchLength=par->GetInnerPadPitchLength();
 932      f2PadPitchLength=f1PadPitchLength;
 933
 934      fN=par->GetNRowLow();
 935      fRow=new AliTPCRow[fN];
 936      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
 937   } else {
 938      fAlpha=par->GetOuterAngle();
 939      fAlphaShift=par->GetOuterAngleShift();
 940      fPadPitchWidth=par->GetOuterPadPitchWidth();
 941      f1PadPitchLength=par->GetOuter1PadPitchLength();
 942      f2PadPitchLength=par->GetOuter2PadPitchLength();
 943
 944      fN=par->GetNRowUp();
 945      fRow=new AliTPCRow[fN];
 946      for (Int_t i=0; i<fN; i++){
 947      fRow[i].SetX(par->GetPadRowRadiiUp(i));
 948      }
 949   }
 950 }
 951
 952 //_________________________________________________________________________
 953 void AliTPCtracker::
 954 AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
 955   //-----------------------------------------------------------------------
 956   // Insert a cluster into this pad row in accordence with its y-coordinate
 957   //-----------------------------------------------------------------------
 958   if (fN==kMaxClusterPerRow) {
 959      ::Error("InsertCluster","Too many clusters !");
 960      return;
 961   }
 962
 963   Int_t index=(((sec<<8)+row)<<16)+fN;
 964
 965   if (fN==0) {
 966      fSize=kMaxClusterPerRow/8;
 967      fClusterArray=new AliTPCcluster[fSize];
 968      fIndex[0]=index;
 969      fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
 970      return;
 971   }
 972
 973   if (fN==fSize) {
 974      Int_t size=fSize*2;
 975      AliTPCcluster *buff=new AliTPCcluster[size];
 976      memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
 977      for (Int_t i=0; i<fN; i++)
 978         fClusters[i]=buff+(fClusters[i]-fClusterArray);
 979      delete[] fClusterArray;
 980      fClusterArray=buff;
 981      fSize=size;
 982   }
 983
 984   Int_t i=Find(c->GetY());
 985   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
 986   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
 987   fIndex[i]=index;
 988   fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
 989 }
 990
 991 //___________________________________________________________________
 992 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
 993   //-----------------------------------------------------------------------
 994   // Return the index of the nearest cluster
 995   //-----------------------------------------------------------------------
 996   if(fN<=0) return 0;
 997   if (y <= fClusters[0]->GetY()) return 0;
 998   if (y > fClusters[fN-1]->GetY()) return fN;
 999   Int_t b=0, e=fN-1, m=(b+e)/2;
1000   for (; b<e; m=(b+e)/2) {
1001     if (y > fClusters[m]->GetY()) b=m+1;
1002     else e=m;
1003   }
1004   return m;
1005 }
1006
1007 //_____________________________________________________________________________
1008 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1009   //-----------------------------------------------------------------
1010   // This funtion calculates dE/dX within the "low" and "up" cuts.
1011   //-----------------------------------------------------------------
1012   Int_t i;
1013   Int_t nc=GetNumberOfClusters();
1014
1015   Int_t swap;//stupid sorting
1016   do {
1017     swap=0;
1018     for (i=0; i<nc-1; i++) {
1019       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1020       Float_t tmp=fdEdxSample[i];
1021       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1022       swap++;
1023     }
1024   } while (swap);
1025
1026   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
1027   Float_t dedx=0;
1028   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
1029   dedx /= (nu-nl+1);
1030   SetdEdx(dedx);
1031
1032   return;
1033 }
1034
1035
1036