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