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 /*
  17 $Log$
  18 Revision 1.30  2003/02/28 16:13:32  hristov
  19 Typos corrected
  20
  21 Revision 1.29  2003/02/28 15:18:16  hristov
  22 Corrections suggested by J.Chudoba
  23
  24 Revision 1.28  2003/02/27 16:15:52  hristov
  25 Code for inward refitting (S.Radomski)
  26
  27 Revision 1.27  2003/02/25 16:47:58  hristov
  28 allow back propagation for more than 1 event (J.Chudoba)
  29
  30 Revision 1.26  2003/02/19 08:49:46  hristov
  31 Track time measurement (S.Radomski)
  32
  33 Revision 1.25  2003/01/28 16:43:35  hristov
  34 Additional protection: to be discussed with the Root team (M.Ivanov)
  35
  36 Revision 1.24  2002/11/19 16:13:24  hristov
  37 stdlib.h included to declare exit() on HP
  38
  39 Revision 1.23  2002/11/19 11:50:08  hristov
  40 Removing CONTAINERS (Yu.Belikov)
  41
  42 Revision 1.19  2002/07/19 07:31:40  kowal2
  43 Improvement in tracking by J. Belikov
  44
  45 Revision 1.18  2002/05/13 07:33:52  kowal2
  46 Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
  47 in the case of defined region of interests.
  48
  49 Revision 1.17  2002/03/18 17:59:13  kowal2
  50 Chnges in the pad geometry - 3 pad lengths introduced.
  51
  52 Revision 1.16  2001/11/08 16:39:03  hristov
  53 Additional protection (M.Masera)
  54
  55 Revision 1.15  2001/11/08 16:36:33  hristov
  56 Updated V2 stream of tracking (Yu.Belikov). The new long waited features are: 1) Possibility to pass the primary vertex position to the trackers (both for the TPC and the ITS) 2) Possibility to specify the number of tracking passes together with applying (or not applying) the vertex constraint (ITS only) 3) Possibility to make some use of partial PID provided by the TPC when doing tracking in the ITS (ITS only) 4) V0 reconstruction with a helix minimisation of the DCA. (new macros: AliV0FindVertices.C and AliV0Comparison.C) 4a) ( Consequence of the 4) )  All the efficiencies and resolutions are from now on calculated including *secondary*particles* too. (Don't be surprised by the drop in efficiency etc)
  57
  58 Revision 1.14  2001/10/21 19:04:55  hristov
  59 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
  60
  61 Revision 1.13  2001/07/23 12:01:30  hristov
  62 Initialisation added
  63
  64 Revision 1.12  2001/07/20 14:32:44  kowal2
  65 Processing of many events possible now
  66
  67 Revision 1.11  2001/05/23 08:50:10  hristov
  68 Weird inline removed
  69
  70 Revision 1.10  2001/05/16 14:57:25  alibrary
  71 New files for folders and Stack
  72
  73 Revision 1.9  2001/05/11 07:16:56  hristov
  74 Fix needed on Sun and Alpha
  75
  76 Revision 1.8  2001/05/08 15:00:15  hristov
  77 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
  78
  79 Revision 1.5  2000/12/20 07:51:59  kowal2
  80 Changes suggested by Alessandra and Paolo to avoid overlapped
  81 data fields in encapsulated classes.
  82
  83 Revision 1.4  2000/11/02 07:27:16  kowal2
  84 code corrections
  85
  86 Revision 1.2  2000/06/30 12:07:50  kowal2
  87 Updated from the TPC-PreRelease branch
  88
  89 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
  90 Splitted from AliTPCtracking
  91
  92 */
  93
  94 //-------------------------------------------------------
  95 //          Implementation of the TPC tracker
  96 //
  97 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  98 //-------------------------------------------------------
  99 #include <TObjArray.h>
 100 #include <TFile.h>
 101 #include <TTree.h>
 102 #include <Riostream.h>
 103
 104 #include "AliTPCtracker.h"
 105 #include "AliTPCcluster.h"
 106 #include "AliTPCParam.h"
 107 #include "AliClusters.h"
 108
 109 #include "TVector2.h"
 110
 111 #include <stdlib.h>
 112 //_____________________________________________________________________________
 113 AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
 114 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 115 {
 116   //---------------------------------------------------------------------
 117   // The main TPC tracker constructor
 118   //---------------------------------------------------------------------
 119   fInnerSec=new AliTPCSector[fkNIS];
 120   fOuterSec=new AliTPCSector[fkNOS];
 121
 122   Int_t i;
 123   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
 124   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
 125
 126   fParam = (AliTPCParam*) par;
 127   fSeeds=0;
 128 }
 129
 130 //_____________________________________________________________________________
 131 AliTPCtracker::~AliTPCtracker() {
 132   //------------------------------------------------------------------
 133   // TPC tracker destructor
 134   //------------------------------------------------------------------
 135   delete[] fInnerSec;
 136   delete[] fOuterSec;
 137   if (fSeeds) {
 138     fSeeds->Delete();
 139     delete fSeeds;
 140   }
 141 }
 142
 143 //_____________________________________________________________________________
 144 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
 145 {
 146   //
 147   // Parametrised error of the cluster reconstruction (pad direction)
 148   //
 149   // Sigma rphi
 150   const Float_t kArphi=0.41818e-2;
 151   const Float_t kBrphi=0.17460e-4;
 152   const Float_t kCrphi=0.30993e-2;
 153   const Float_t kDrphi=0.41061e-3;
 154
 155   pt=TMath::Abs(pt)*1000.;
 156   Double_t x=r/pt;
 157   tgl=TMath::Abs(tgl);
 158   Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
 159   if (s<0.4e-3) s=0.4e-3;
 160   s*=1.3; //Iouri Belikov
 161
 162   return s;
 163 }
 164
 165 //_____________________________________________________________________________
 166 Double_t SigmaZ2(Double_t r, Double_t tgl)
 167 {
 168   //
 169   // Parametrised error of the cluster reconstruction (drift direction)
 170   //
 171   // Sigma z
 172   const Float_t kAz=0.39614e-2;
 173   const Float_t kBz=0.22443e-4;
 174   const Float_t kCz=0.51504e-1;
 175
 176
 177   tgl=TMath::Abs(tgl);
 178   Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
 179   if (s<0.4e-3) s=0.4e-3;
 180   s*=1.3; //Iouri Belikov
 181
 182   return s;
 183 }
 184
 185 //_____________________________________________________________________________
 186 Double_t f1(Double_t x1,Double_t y1,
 187                    Double_t x2,Double_t y2,
 188                    Double_t x3,Double_t y3)
 189 {
 190   //-----------------------------------------------------------------
 191   // Initial approximation of the track curvature
 192   //-----------------------------------------------------------------
 193   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 194   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 195                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 196   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 197                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 198
 199   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 200
 201   return -xr*yr/sqrt(xr*xr+yr*yr);
 202 }
 203
 204
 205 //_____________________________________________________________________________
 206 Double_t f2(Double_t x1,Double_t y1,
 207                    Double_t x2,Double_t y2,
 208                    Double_t x3,Double_t y3)
 209 {
 210   //-----------------------------------------------------------------
 211   // Initial approximation of the track curvature times center of curvature
 212   //-----------------------------------------------------------------
 213   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 214   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 215                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 216   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 217                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 218
 219   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 220
 221   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 222 }
 223
 224 //_____________________________________________________________________________
 225 Double_t f3(Double_t x1,Double_t y1,
 226                    Double_t x2,Double_t y2,
 227                    Double_t z1,Double_t z2)
 228 {
 229   //-----------------------------------------------------------------
 230   // Initial approximation of the tangent of the track dip angle
 231   //-----------------------------------------------------------------
 232   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 233 }
 234
 235 //_____________________________________________________________________________
 236 Int_t AliTPCtracker::LoadClusters() {
 237   //-----------------------------------------------------------------
 238   // This function loads TPC clusters.
 239   //-----------------------------------------------------------------
 240   if (!gFile->IsOpen()) {
 241     cerr<<"AliTPCtracker::LoadClusters : "<<
 242       "file with clusters has not been open !\n";
 243     return 1;
 244   }
 245
 246   Char_t name[100];
 247   sprintf(name,"TreeC_TPC_%d",GetEventNumber());
 248   TTree *cTree=(TTree*)gFile->Get(name);
 249   if (!cTree) {
 250     cerr<<"AliTPCtracker::LoadClusters : "<<
 251       "can't get the tree with TPC clusters !\n";
 252     return 2;
 253   }
 254
 255   TBranch *branch=cTree->GetBranch("Segment");
 256   if (!branch) {
 257     cerr<<"AliTPCtracker::LoadClusters : "<<
 258       "can't get the segment branch !\n";
 259     return 3;
 260   }
 261 //  AliClusters carray, *addr=&carray;
 262   AliClusters carray, *addr=&carray;
 263   carray.SetClass("AliTPCcluster");
 264   carray.SetArray(0);
 265   branch->SetAddress(&addr);
 266
 267   Int_t nentr=(Int_t)cTree->GetEntries();
 268
 269   for (Int_t i=0; i<nentr; i++) {
 270       cTree->GetEvent(i);
 271
 272       Int_t ncl=carray.GetArray()->GetEntriesFast();
 273
 274       Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
 275       Int_t id=carray.GetID();
 276       if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
 277          cerr<<"AliTPCtracker::LoadClusters : "<<
 278                "wrong index !\n";
 279          exit(1);
 280       }
 281       Int_t outindex = 2*fkNIS*nir;
 282       if (id<outindex) {
 283          Int_t sec = id/nir;
 284          Int_t row = id - sec*nir;
 285          sec %= fkNIS;
 286          AliTPCRow &padrow=fInnerSec[sec][row];
 287          while (ncl--) {
 288            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
 289            padrow.InsertCluster(c,sec,row);
 290          }
 291       } else {
 292          id -= outindex;
 293          Int_t sec = id/nor;
 294          Int_t row = id - sec*nor;
 295          sec %= fkNOS;
 296          AliTPCRow &padrow=fOuterSec[sec][row];
 297          while (ncl--) {
 298            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
 299            padrow.InsertCluster(c,sec+fkNIS,row);
 300          }
 301       }
 302
 303       carray.GetArray()->Clear();
 304   }
 305   delete cTree;
 306   return 0;
 307 }
 308
 309 //_____________________________________________________________________________
 310 void AliTPCtracker::UnloadClusters() {
 311   //-----------------------------------------------------------------
 312   // This function unloads TPC clusters.
 313   //-----------------------------------------------------------------
 314   Int_t i;
 315   for (i=0; i<fkNIS; i++) {
 316     Int_t nr=fInnerSec->GetNRows();
 317     for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
 318   }
 319   for (i=0; i<fkNOS; i++) {
 320     Int_t nr=fOuterSec->GetNRows();
 321     for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
 322   }
 323 }
 324
 325 //_____________________________________________________________________________
 326 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
 327   //-----------------------------------------------------------------
 328   // This function tries to find a track prolongation.
 329   //-----------------------------------------------------------------
 330   Double_t xt=t.GetX();
 331   const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
 332                     Int_t(0.5*fSectors->GetNRows());
 333   Int_t tryAgain=kSKIP;
 334
 335   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
 336   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 337   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 338   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 339
 340   Int_t nrows=fSectors->GetRowNumber(xt)-1;
 341   for (Int_t nr=nrows; nr>=rf; nr--) {
 342     Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
 343     if (!t.PropagateTo(x)) return 0;
 344
 345     AliTPCcluster *cl=0;
 346     UInt_t index=0;
 347     Double_t maxchi2=kMaxCHI2;
 348     const AliTPCRow &krow=fSectors[s][nr];
 349     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
 350     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
 351     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
 352     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 353
 354     if (road>kMaxROAD) {
 355       if (t.GetNumberOfClusters()>4)
 356         cerr<<t.GetNumberOfClusters()
 357         <<"FindProlongation warning: Too broad road !\n";
 358       return 0;
 359     }
 360
 361     if (krow) {
 362       for (Int_t i=krow.Find(y-road); i<krow; i++) {
 363         AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 364         if (c->GetY() > y+road) break;
 365         if (c->IsUsed()) continue;
 366         if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
 367         Double_t chi2=t.GetPredictedChi2(c);
 368         if (chi2 > maxchi2) continue;
 369         maxchi2=chi2;
 370         cl=c;
 371         index=krow.GetIndex(i);
 372       }
 373     }
 374     if (cl) {
 375       Float_t l=fSectors->GetPadPitchWidth();
 376       Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
 377       t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
 378       if (!t.Update(cl,maxchi2,index)) {
 379          if (!tryAgain--) return 0;
 380       } else tryAgain=kSKIP;
 381     } else {
 382       if (tryAgain==0) break;
 383       if (y > ymax) {
 384          s = (s+1) % fN;
 385          if (!t.Rotate(fSectors->GetAlpha())) return 0;
 386       } else if (y <-ymax) {
 387          s = (s-1+fN) % fN;
 388          if (!t.Rotate(-fSectors->GetAlpha())) return 0;
 389       }
 390       tryAgain--;
 391     }
 392   }
 393
 394   return 1;
 395 }
 396 //_____________________________________________________________________________
 397
 398 Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
 399   //
 400   // This function propagates seed inward TPC using old clusters
 401   // from track.
 402   //
 403   // Sylwester Radomski, GSI
 404   // 26.02.2003
 405   //
 406
 407   Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
 408   TVector2::Phi_0_2pi(alpha);
 409   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 410
 411   Int_t idx=-1, sec=-1, row=-1;
 412   Int_t nc = seed->GetLabel(); // index to start with
 413
 414   idx=track->GetClusterIndex(nc);
 415   sec=(idx&0xff000000)>>24;
 416   row=(idx&0x00ff0000)>>16;
 417
 418   if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 419   else { if (sec <  fkNIS) row=-1; }
 420
 421   Int_t nr=fSectors->GetNRows();
 422   for (Int_t i=0; i<nr; i++) {
 423
 424     Double_t x=fSectors->GetX(i);
 425     if (!seed->PropagateTo(x)) return 0;
 426
 427     // Change sector and rotate seed if necessary
 428
 429     Double_t ymax=fSectors->GetMaxY(i);
 430     Double_t y=seed->GetY();
 431
 432     if (y > ymax) {
 433        s = (s+1) % fN;
 434        if (!seed->Rotate(fSectors->GetAlpha())) return 0;
 435     } else if (y <-ymax) {
 436        s = (s-1+fN) % fN;
 437        if (!seed->Rotate(-fSectors->GetAlpha())) return 0;
 438     }
 439
 440     // try to find a cluster
 441
 442     AliTPCcluster *cl=0;
 443     Int_t index = 0;
 444     Double_t maxchi2 = kMaxCHI2;
 445
 446     //cout << i << " " << row << " " << nc << endl;
 447
 448     if (row==i) {
 449
 450       // accept already found cluster
 451       AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
 452       Double_t chi2 = seed->GetPredictedChi2(c);
 453
 454       index=idx;
 455       cl=c;
 456       maxchi2=chi2;
 457
 458       if (++nc < track->GetNumberOfClusters() ) {
 459         idx=track->GetClusterIndex(nc);
 460         sec=(idx&0xff000000)>>24;
 461         row=(idx&0x00ff0000)>>16;
 462       }
 463
 464       if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 465       else { if (sec < fkNIS) row=-1; }
 466
 467     }
 468
 469     if (cl) {
 470
 471       //cout << "Assigned" << endl;
 472       Float_t l=fSectors->GetPadPitchWidth();
 473       Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
 474       seed->SetSampledEdx(cl->GetQ()/l*corr,seed->GetNumberOfClusters());
 475       seed->Update(cl,maxchi2,index);
 476     }
 477   }
 478
 479   seed->SetLabel(nc);
 480   return 1;
 481 }
 482
 483 //_____________________________________________________________________________
 484 Int_t AliTPCtracker::FollowBackProlongation
 485 (AliTPCseed& seed, const AliTPCtrack &track) {
 486   //-----------------------------------------------------------------
 487   // This function propagates tracks back through the TPC
 488   //-----------------------------------------------------------------
 489   Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
 490   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 491   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 492   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 493
 494   Int_t idx=-1, sec=-1, row=-1;
 495   Int_t nc=seed.GetLabel(); //index of the cluster to start with
 496   if (nc--) {
 497      idx=track.GetClusterIndex(nc);
 498      sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 499   }
 500   if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 501   else { if (sec <  fkNIS) row=-1; }
 502
 503   Int_t nr=fSectors->GetNRows();
 504   for (Int_t i=0; i<nr; i++) {
 505     Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
 506
 507     if (!seed.PropagateTo(x)) return 0;
 508
 509     Double_t y=seed.GetY();
 510     if (y > ymax) {
 511        s = (s+1) % fN;
 512        if (!seed.Rotate(fSectors->GetAlpha())) return 0;
 513     } else if (y <-ymax) {
 514        s = (s-1+fN) % fN;
 515        if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
 516     }
 517
 518     AliTPCcluster *cl=0;
 519     Int_t index=0;
 520     Double_t maxchi2=kMaxCHI2;
 521     Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
 522     Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
 523     Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
 524     Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
 525     if (road>kMaxROAD) {
 526       cerr<<seed.GetNumberOfClusters()
 527           <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
 528       return 0;
 529     }
 530
 531
 532     Int_t accepted=seed.GetNumberOfClusters();
 533     if (row==i) {
 534        //try to accept already found cluster
 535        AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
 536        Double_t chi2;
 537        if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
 538           index=idx; cl=c; maxchi2=chi2;
 539        } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
 540
 541        if (nc--) {
 542           idx=track.GetClusterIndex(nc);
 543           sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 544        }
 545        if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 546        else { if (sec < fkNIS) row=-1; }
 547
 548     }
 549     if (!cl) {
 550        //try to fill the gap
 551        const AliTPCRow &krow=fSectors[s][i];
 552        if (accepted>27)
 553        if (krow) {
 554           for (Int_t i=krow.Find(y-road); i<krow; i++) {
 555             AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 556             if (c->GetY() > y+road) break;
 557             if (c->IsUsed()) continue;
 558          if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
 559             Double_t chi2=seed.GetPredictedChi2(c);
 560             if (chi2 > maxchi2) continue;
 561             maxchi2=chi2;
 562             cl=c;
 563             index=krow.GetIndex(i);
 564           }
 565        }
 566     }
 567
 568     if (cl) {
 569       Float_t l=fSectors->GetPadPitchWidth();
 570       Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
 571       seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
 572       seed.Update(cl,maxchi2,index);
 573     }
 574
 575   }
 576
 577   seed.SetLabel(nc);
 578
 579   return 1;
 580 }
 581
 582 //_____________________________________________________________________________
 583 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
 584   //-----------------------------------------------------------------
 585   // This function creates track seeds.
 586   //-----------------------------------------------------------------
 587   Double_t x[5], c[15];
 588
 589   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
 590   Double_t cs=cos(alpha), sn=sin(alpha);
 591
 592   Double_t x1 =fSectors->GetX(i1);
 593   Double_t xx2=fSectors->GetX(i2);
 594
 595   for (Int_t ns=0; ns<fN; ns++) {
 596     Int_t nl=fSectors[(ns-1+fN)%fN][i2];
 597     Int_t nm=fSectors[ns][i2];
 598     Int_t nu=fSectors[(ns+1)%fN][i2];
 599     const AliTPCRow& kr1=fSectors[ns][i1];
 600     for (Int_t is=0; is < kr1; is++) {
 601       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
 602       for (Int_t js=0; js < nl+nm+nu; js++) {
 603         const AliTPCcluster *kcl;
 604         Double_t x2,   y2,   z2;
 605         Double_t x3=GetX(), y3=GetY(), z3=GetZ();
 606
 607         if (js<nl) {
 608           const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
 609           kcl=kr2[js];
 610           y2=kcl->GetY(); z2=kcl->GetZ();
 611           x2= xx2*cs+y2*sn;
 612           y2=-xx2*sn+y2*cs;
 613         } else
 614           if (js<nl+nm) {
 615             const AliTPCRow& kr2=fSectors[ns][i2];
 616             kcl=kr2[js-nl];
 617             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
 618           } else {
 619             const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
 620             kcl=kr2[js-nl-nm];
 621             y2=kcl->GetY(); z2=kcl->GetZ();
 622             x2=xx2*cs-y2*sn;
 623             y2=xx2*sn+y2*cs;
 624           }
 625
 626         Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
 627         if (TMath::Abs(zz-z2)>5.) continue;
 628
 629         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
 630         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
 631
 632         x[0]=y1;
 633         x[1]=z1;
 634         x[4]=f1(x1,y1,x2,y2,x3,y3);
 635         if (TMath::Abs(x[4]) >= 0.0066) continue;
 636         x[2]=f2(x1,y1,x2,y2,x3,y3);
 637         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
 638         x[3]=f3(x1,y1,x2,y2,z1,z2);
 639         if (TMath::Abs(x[3]) > 1.2) continue;
 640         Double_t a=asin(x[2]);
 641         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
 642         if (TMath::Abs(zv-z3)>10.) continue;
 643
 644         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
 645         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
 646         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
 647         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
 648
 649         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
 650         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
 651         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
 652         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
 653         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
 654         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
 655         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
 656         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
 657         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
 658         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
 659
 660         c[0]=sy1;
 661         c[1]=0.;       c[2]=sz1;
 662         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
 663         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
 664                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
 665         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
 666         c[13]=f30*sy1*f40+f32*sy2*f42;
 667         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
 668
 669         UInt_t index=kr1.GetIndex(is);
 670         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
 671         Float_t l=fSectors->GetPadPitchWidth();
 672         track->SetSampledEdx(kr1[is]->GetQ()/l,0);
 673
 674         Int_t rc=FollowProlongation(*track, i2);
 675         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
 676         else fSeeds->AddLast(track);
 677       }
 678     }
 679   }
 680 }
 681
 682 //_____________________________________________________________________________
 683 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
 684   //-----------------------------------------------------------------
 685   // This function reades track seeds.
 686   //-----------------------------------------------------------------
 687   TDirectory *savedir=gDirectory;
 688
 689   TFile *in=(TFile*)inp;
 690   if (!in->IsOpen()) {
 691      cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
 692      return 1;
 693   }
 694
 695   in->cd();
 696   TTree *seedTree=(TTree*)in->Get("Seeds");
 697   if (!seedTree) {
 698      cerr<<"AliTPCtracker::ReadSeeds(): ";
 699      cerr<<"can't get a tree with track seeds !\n";
 700      return 2;
 701   }
 702   AliTPCtrack *seed=new AliTPCtrack;
 703   seedTree->SetBranchAddress("tracks",&seed);
 704
 705   if (fSeeds==0) fSeeds=new TObjArray(15000);
 706
 707   Int_t n=(Int_t)seedTree->GetEntries();
 708   for (Int_t i=0; i<n; i++) {
 709      seedTree->GetEvent(i);
 710      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
 711   }
 712
 713   delete seed;
 714
 715   delete seedTree; //Thanks to Mariana Bondila
 716
 717   savedir->cd();
 718
 719   return 0;
 720 }
 721
 722 //_____________________________________________________________________________
 723 Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
 724   //-----------------------------------------------------------------
 725   // This is a track finder.
 726   //-----------------------------------------------------------------
 727   TDirectory *savedir=gDirectory;
 728
 729   if (inp) {
 730      TFile *in=(TFile*)inp;
 731      if (!in->IsOpen()) {
 732         cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
 733         return 1;
 734      }
 735   }
 736
 737   if (!out->IsOpen()) {
 738      cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
 739      return 2;
 740   }
 741
 742   LoadClusters();
 743
 744   out->cd();
 745
 746   char   tname[100];
 747   sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
 748   TTree tracktree(tname,"Tree with TPC tracks");
 749   AliTPCtrack *iotrack=0;
 750   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,3);
 751
 752   //find track seeds
 753   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 754   Int_t nrows=nlow+nup;
 755   if (fSeeds==0) {
 756      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
 757      fSectors=fOuterSec; fN=fkNOS;
 758      fSeeds=new TObjArray(15000);
 759      MakeSeeds(nup-1, nup-1-gap);
 760      MakeSeeds(nup-1-shift, nup-1-shift-gap);
 761   }
 762   fSeeds->Sort();
 763
 764   Int_t found=0;
 765   Int_t nseed=fSeeds->GetEntriesFast();
 766   for (Int_t i=0; i<nseed; i++) {
 767     //tracking in the outer sectors
 768     fSectors=fOuterSec; fN=fkNOS;
 769
 770     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 771     if (!FollowProlongation(t)) {
 772        delete fSeeds->RemoveAt(i);
 773        continue;
 774     }
 775
 776     //tracking in the inner sectors
 777     fSectors=fInnerSec; fN=fkNIS;
 778
 779     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
 780     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 781     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 782     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
 783
 784     alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
 785
 786     if (t.Rotate(alpha)) {
 787       if (FollowProlongation(t)) {
 788         if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 789           t.CookdEdx();
 790           CookLabel(pt,0.1); //For comparison only
 791           pt->PropagateTo(fParam->GetInnerRadiusLow());
 792           iotrack=pt;
 793           tracktree.Fill();
 794           UseClusters(&t);
 795           found++;
 796 //          cerr<<found<<'\r';
 797         }
 798       }
 799     }
 800     delete fSeeds->RemoveAt(i);
 801   }
 802
 803   tracktree.Write();
 804
 805   cerr<<"Number of found tracks : "<<found<<endl;
 806
 807   savedir->cd();
 808
 809   UnloadClusters();
 810   fSeeds->Clear(); delete fSeeds; fSeeds=0;
 811
 812   return 0;
 813 }
 814 //_____________________________________________________________________________
 815
 816 Int_t AliTPCtracker::RefitInward(TFile *in, TFile *out) {
 817   //
 818   // The function propagates tracks throught TPC inward
 819   // using already associated clusters.
 820   //
 821   // in - file with back propagated tracks
 822   // outs  - file with inward propagation
 823   //
 824   // Sylwester Radomski, GSI
 825   // 26.02.2003
 826   //
 827
 828
 829   if (!in->IsOpen()) {
 830     cout << "Input File not open !\n" << endl;
 831     return 1;
 832   }
 833
 834   if (!out->IsOpen()) {
 835     cout << "Output File not open !\n" << endl;
 836     return 2;
 837   }
 838
 839   TDirectory *savedir = gDirectory;
 840   LoadClusters();
 841
 842   // Connect to input tree
 843   in->cd();
 844   TTree *inputTree = (TTree*)in->Get("seedsTPCtoTRD_0");
 845   TBranch *inBranch = inputTree->GetBranch("tracks");
 846   AliTPCtrack *inTrack = 0;
 847   inBranch->SetAddress(&inTrack);
 848
 849   Int_t nTracks = (Int_t)inputTree->GetEntries();
 850
 851   // Create output tree
 852   out->cd();
 853   AliTPCtrack *outTrack = new AliTPCtrack();
 854   TTree *outputTree = new TTree("tracksTPC_0","Refited TPC tracks");
 855   outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
 856
 857   Int_t nRefited = 0;
 858
 859   for(Int_t t=0; t < nTracks; t++) {
 860
 861     cout << t << "\r";
 862
 863     inputTree->GetEvent(t);
 864     AliTPCseed *seed = new AliTPCseed(*inTrack, inTrack->GetAlpha());
 865
 866     seed->ResetCovariance();
 867
 868     seed->SetLabel(0);
 869     fSectors = fInnerSec;
 870     Int_t res = FollowRefitInward(seed, inTrack);
 871     UseClusters(seed);
 872     Int_t nc = seed->GetNumberOfClusters();
 873
 874     fSectors = fOuterSec;
 875     res = FollowRefitInward(seed, inTrack);
 876     UseClusters(seed, nc);
 877
 878     if (res) {
 879       seed->PropagateTo(fParam->GetInnerRadiusLow());
 880       outTrack = (AliTPCtrack*)seed;
 881       outTrack->SetLabel(inTrack->GetLabel());
 882       outputTree->Fill();
 883       nRefited++;
 884     }
 885
 886     delete seed;
 887   }
 888
 889   cout << "Refitted " << nRefited << " tracks." << endl;
 890
 891   outputTree->Write();
 892
 893   if (inputTree) delete inputTree;
 894   if (outputTree) delete outputTree;
 895
 896   savedir->cd();
 897   UnloadClusters();
 898
 899   return 0;
 900 }
 901
 902 //_____________________________________________________________________________
 903 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
 904   //-----------------------------------------------------------------
 905   // This function propagates tracks back through the TPC.
 906   //-----------------------------------------------------------------
 907   return PropagateBack(inp, NULL, out);
 908 }
 909
 910 //_____________________________________________________________________________
 911 Int_t AliTPCtracker::PropagateBack(const TFile *inp, const TFile *inp2, TFile *out) {
 912   //-----------------------------------------------------------------
 913   // This function propagates tracks back through the TPC.
 914   //-----------------------------------------------------------------
 915   fSeeds=new TObjArray(15000);
 916
 917   TFile *in=(TFile*)inp;
 918   TFile *in2=(TFile*)inp2;
 919   TDirectory *savedir=gDirectory;
 920
 921   if (!in->IsOpen()) {
 922      cerr<<"AliTPCtracker::PropagateBack(): ";
 923      cerr<<"file with TPC (or back propagated ITS) tracks is not open !\n";
 924      return 1;
 925   }
 926
 927   if (!out->IsOpen()) {
 928      cerr<<"AliTPCtracker::PropagateBack(): ";
 929      cerr<<"file for back propagated TPC tracks is not open !\n";
 930      return 2;
 931   }
 932
 933   LoadClusters();
 934
 935   in->cd();
 936   char tName[100];
 937   sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
 938   TTree *bckTree=(TTree*)in->Get(tName);
 939   if (!bckTree && inp2) bckTree=(TTree*)in2->Get(tName);
 940   if (!bckTree) {
 941      cerr<<"AliTPCtracker::PropagateBack() ";
 942      cerr<<"can't get a tree with back propagated ITS tracks !\n";
 943      //return 3;
 944   }
 945   AliTPCtrack *bckTrack=new AliTPCtrack;
 946   if (bckTree) bckTree->SetBranchAddress("tracks",&bckTrack);
 947
 948   sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
 949   TTree *tpcTree=(TTree*)in->Get(tName);
 950   if (!tpcTree) {
 951      cerr<<"AliTPCtracker::PropagateBack() ";
 952      cerr<<"can't get a tree with TPC tracks !\n";
 953      return 4;
 954   }
 955   AliTPCtrack *tpcTrack=new AliTPCtrack;
 956   tpcTree->SetBranchAddress("tracks",&tpcTrack);
 957
 958   //*** Prepare an array of tracks to be back propagated
 959   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 960   Int_t nrows=nlow+nup;
 961
 962   //
 963   // Match ITS tracks with old TPC tracks
 964   // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
 965   //
 966   // the algorithm is linear and uses LUT for sorting
 967   // cost of the algorithm: 2 * number of tracks
 968   //
 969
 970   TObjArray tracks(15000);
 971   Int_t i=0;
 972   Int_t tpcN= (Int_t)tpcTree->GetEntries();
 973   Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
 974
 975   // look up table
 976   const Int_t nLab = 200000; // limit on number of primaries (arbitrary)
 977   Int_t lut[nLab];
 978   for(Int_t i=0; i<nLab; i++) lut[i] = -1;
 979
 980   if (bckTree) {
 981     for(Int_t i=0; i<bckN; i++) {
 982       bckTree->GetEvent(i);
 983       Int_t lab = TMath::Abs(bckTrack->GetLabel());
 984       if (lab < nLab) lut[lab] = i;
 985       else {
 986         cerr << "AliTPCtracker: The size of the LUT is too small\n";
 987       }
 988     }
 989   }
 990
 991   for (Int_t i=0; i<tpcN; i++) {
 992     tpcTree->GetEvent(i);
 993     Double_t alpha=tpcTrack->GetAlpha();
 994
 995     if (!bckTree) {
 996
 997       // No ITS - use TPC track only
 998
 999       AliTPCseed * seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
1000       seed->ResetCovariance();
1001       fSeeds->AddLast(seed);
1002       tracks.AddLast(new AliTPCtrack(*tpcTrack));
1003
1004     } else {
1005
1006       // with ITS
1007       // discard not prolongated tracks (to be discussed)
1008
1009       Int_t lab = TMath::Abs(tpcTrack->GetLabel());
1010       if (lab > nLab || lut[lab] < 0) continue;
1011
1012       bckTree->GetEvent(lut[lab]);
1013       bckTrack->Rotate(alpha-bckTrack->GetAlpha());
1014
1015       fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1016       tracks.AddLast(new AliTPCtrack(*tpcTrack));
1017     }
1018   }
1019
1020   // end of matching
1021
1022   /*
1023   TObjArray tracks(15000);
1024   Int_t i=0, j=0;
1025   Int_t tpcN=(Int_t)tpcTree->GetEntries();
1026   Int_t bckN=(Int_t)bckTree->GetEntries();
1027   if (j<bckN) bckTree->GetEvent(j++);
1028   for (i=0; i<tpcN; i++) {
1029      tpcTree->GetEvent(i);
1030      Double_t alpha=tpcTrack->GetAlpha();
1031      if (j<bckN &&
1032      TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
1033         if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
1034         fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1035         bckTree->GetEvent(j++);
1036      } else {
1037         tpcTrack->ResetCovariance();
1038         fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
1039      }
1040      tracks.AddLast(new AliTPCtrack(*tpcTrack));
1041   }
1042   */
1043
1044   out->cd();
1045
1046   // tree name seedsTPCtoTRD as expected by TRD and as
1047   // discussed and decided in Strasbourg (May 2002)
1048   // [SR, GSI, 18.02.2003]
1049
1050   sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
1051   TTree backTree(tName,"Tree with back propagated TPC tracks");
1052
1053   AliTPCtrack *otrack=0;
1054   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
1055
1056   //*** Back propagation through inner sectors
1057   fSectors=fInnerSec; fN=fkNIS;
1058   Int_t nseed=fSeeds->GetEntriesFast();
1059   for (i=0; i<nseed; i++) {
1060      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1061      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1062
1063      Int_t nc=t.GetNumberOfClusters();
1064      s.SetLabel(nc-1); //set number of the cluster to start with
1065
1066      if (FollowBackProlongation(s,t)) {
1067         UseClusters(&s);
1068         continue;
1069      }
1070      delete fSeeds->RemoveAt(i);
1071   }
1072
1073 //*** Back propagation through outer sectors
1074   fSectors=fOuterSec; fN=fkNOS;
1075   Int_t found=0;
1076   for (i=0; i<nseed; i++) {
1077      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1078      if (!ps) continue;
1079      Int_t nc=s.GetNumberOfClusters();
1080      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1081
1082      Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1083      if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1084      if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1085      Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1086
1087      alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1088      alpha-=s.GetAlpha();
1089
1090      if (s.Rotate(alpha)) {
1091        if (FollowBackProlongation(s,t)) {
1092          if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
1093            s.CookdEdx();
1094            s.SetLabel(t.GetLabel());
1095            UseClusters(&s,nc);
1096
1097            // Propagate to outer reference plane for comparison reasons
1098            // reason for keeping fParam object [SR, GSI, 18.02.2003]
1099
1100            ps->PropagateTo(fParam->GetOuterRadiusUp());
1101            otrack=ps;
1102            backTree.Fill();
1103            found++;
1104 //           cerr<<found<<'\r';
1105            continue;
1106          }
1107        }
1108      }
1109      delete fSeeds->RemoveAt(i);
1110   }
1111
1112   backTree.Write();
1113   savedir->cd();
1114   cerr<<"Number of seeds: "<<nseed<<endl;
1115   cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
1116   cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
1117
1118   delete bckTrack;
1119   delete tpcTrack;
1120
1121   if (bckTree) delete bckTree; //Thanks to Mariana Bondila
1122   delete tpcTree; //Thanks to Mariana Bondila
1123
1124   UnloadClusters();
1125
1126   return 0;
1127 }
1128
1129 //_________________________________________________________________________
1130 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
1131   //--------------------------------------------------------------------
1132   //       Return pointer to a given cluster
1133   //--------------------------------------------------------------------
1134   Int_t sec=(index&0xff000000)>>24;
1135   Int_t row=(index&0x00ff0000)>>16;
1136   Int_t ncl=(index&0x0000ffff)>>00;
1137
1138   const AliTPCcluster *cl=0;
1139
1140   if (sec<fkNIS) {
1141     cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
1142   } else {
1143     sec-=fkNIS;
1144     cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
1145   }
1146
1147   return (AliCluster*)cl;
1148 }
1149
1150 //__________________________________________________________________________
1151 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1152   //--------------------------------------------------------------------
1153   //This function "cooks" a track label. If label<0, this track is fake.
1154   //--------------------------------------------------------------------
1155   Int_t noc=t->GetNumberOfClusters();
1156   Int_t *lb=new Int_t[noc];
1157   Int_t *mx=new Int_t[noc];
1158   AliCluster **clusters=new AliCluster*[noc];
1159
1160   Int_t i;
1161   for (i=0; i<noc; i++) {
1162      lb[i]=mx[i]=0;
1163      Int_t index=t->GetClusterIndex(i);
1164      clusters[i]=GetCluster(index);
1165   }
1166
1167   Int_t lab=123456789;
1168   for (i=0; i<noc; i++) {
1169     AliCluster *c=clusters[i];
1170     lab=TMath::Abs(c->GetLabel(0));
1171     Int_t j;
1172     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1173     lb[j]=lab;
1174     (mx[j])++;
1175   }
1176
1177   Int_t max=0;
1178   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1179
1180   for (i=0; i<noc; i++) {
1181     AliCluster *c=clusters[i];
1182     if (TMath::Abs(c->GetLabel(1)) == lab ||
1183         TMath::Abs(c->GetLabel(2)) == lab ) max++;
1184   }
1185
1186   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1187
1188   else {
1189      Int_t tail=Int_t(0.10*noc);
1190      max=0;
1191      for (i=1; i<=tail; i++) {
1192        AliCluster *c=clusters[noc-i];
1193        if (lab == TMath::Abs(c->GetLabel(0)) ||
1194            lab == TMath::Abs(c->GetLabel(1)) ||
1195            lab == TMath::Abs(c->GetLabel(2))) max++;
1196      }
1197      if (max < Int_t(0.5*tail)) lab=-lab;
1198   }
1199
1200   t->SetLabel(lab);
1201
1202   delete[] lb;
1203   delete[] mx;
1204   delete[] clusters;
1205 }
1206
1207 //_________________________________________________________________________
1208 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1209   //-----------------------------------------------------------------------
1210   // Setup inner sector
1211   //-----------------------------------------------------------------------
1212   if (f==0) {
1213      fAlpha=par->GetInnerAngle();
1214      fAlphaShift=par->GetInnerAngleShift();
1215      fPadPitchWidth=par->GetInnerPadPitchWidth();
1216      f1PadPitchLength=par->GetInnerPadPitchLength();
1217      f2PadPitchLength=f1PadPitchLength;
1218
1219      fN=par->GetNRowLow();
1220      fRow=new AliTPCRow[fN];
1221      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1222   } else {
1223      fAlpha=par->GetOuterAngle();
1224      fAlphaShift=par->GetOuterAngleShift();
1225      fPadPitchWidth=par->GetOuterPadPitchWidth();
1226      f1PadPitchLength=par->GetOuter1PadPitchLength();
1227      f2PadPitchLength=par->GetOuter2PadPitchLength();
1228
1229      fN=par->GetNRowUp();
1230      fRow=new AliTPCRow[fN];
1231      for (Int_t i=0; i<fN; i++){
1232      fRow[i].SetX(par->GetPadRowRadiiUp(i));
1233      }
1234   }
1235 }
1236
1237 //_________________________________________________________________________
1238 void AliTPCtracker::
1239 AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
1240   //-----------------------------------------------------------------------
1241   // Insert a cluster into this pad row in accordence with its y-coordinate
1242   //-----------------------------------------------------------------------
1243   if (fN==kMaxClusterPerRow) {
1244     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1245   }
1246
1247   Int_t index=(((sec<<8)+row)<<16)+fN;
1248
1249   if (fN==0) {
1250      fIndex[0]=index;
1251      fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1252      return;
1253   }
1254
1255   if (fN==fSize) {
1256      Int_t size=fSize*2;
1257      AliTPCcluster *buff=new AliTPCcluster[size];
1258      memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1259      for (Int_t i=0; i<fN; i++)
1260         fClusters[i]=buff+(fClusters[i]-fClusterArray);
1261      delete[] fClusterArray;
1262      fClusterArray=buff;
1263      fSize=size;
1264   }
1265
1266   Int_t i=Find(c->GetY());
1267   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1268   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
1269   fIndex[i]=index;
1270   fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
1271 }
1272
1273 //___________________________________________________________________
1274 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1275   //-----------------------------------------------------------------------
1276   // Return the index of the nearest cluster
1277   //-----------------------------------------------------------------------
1278   if(fN<=0) return 0;
1279   if (y <= fClusters[0]->GetY()) return 0;
1280   if (y > fClusters[fN-1]->GetY()) return fN;
1281   Int_t b=0, e=fN-1, m=(b+e)/2;
1282   for (; b<e; m=(b+e)/2) {
1283     if (y > fClusters[m]->GetY()) b=m+1;
1284     else e=m;
1285   }
1286   return m;
1287 }
1288
1289 //_____________________________________________________________________________
1290 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1291   //-----------------------------------------------------------------
1292   // This funtion calculates dE/dX within the "low" and "up" cuts.
1293   //-----------------------------------------------------------------
1294   Int_t i;
1295   Int_t nc=GetNumberOfClusters();
1296
1297   Int_t swap;//stupid sorting
1298   do {
1299     swap=0;
1300     for (i=0; i<nc-1; i++) {
1301       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1302       Float_t tmp=fdEdxSample[i];
1303       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1304       swap++;
1305     }
1306   } while (swap);
1307
1308   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
1309   Float_t dedx=0;
1310   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
1311   dedx /= (nu-nl+1);
1312   SetdEdx(dedx);
1313
1314   //Very rough PID
1315   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1316
1317   Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
1318   if (p<0.6) {
1319     if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1320     if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
1321     SetMass(0.93827); return;
1322   }
1323
1324   if (p<1.2) {
1325     if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
1326     SetMass(0.93827); return;
1327   }
1328
1329   SetMass(0.13957); return;
1330
1331 }
1332
1333
1334