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.8  2001/05/08 15:00:15  hristov
  19 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
  20
  21 Revision 1.5  2000/12/20 07:51:59  kowal2
  22 Changes suggested by Alessandra and Paolo to avoid overlapped
  23 data fields in encapsulated classes.
  24
  25 Revision 1.4  2000/11/02 07:27:16  kowal2
  26 code corrections
  27
  28 Revision 1.2  2000/06/30 12:07:50  kowal2
  29 Updated from the TPC-PreRelease branch
  30
  31 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
  32 Splitted from AliTPCtracking
  33
  34 */
  35
  36 //-------------------------------------------------------
  37 //          Implementation of the TPC tracker
  38 //
  39 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  40 //-------------------------------------------------------
  41
  42 #include <TObjArray.h>
  43 #include <TFile.h>
  44 #include <TTree.h>
  45 #include <iostream.h>
  46
  47 #include "AliTPCtracker.h"
  48 #include "AliTPCcluster.h"
  49 #include "AliTPCParam.h"
  50 #include "AliTPCClustersRow.h"
  51
  52 //_____________________________________________________________________________
  53 AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
  54 fkNIS(par->GetNInnerSector()/2),
  55 fkNOS(par->GetNOuterSector()/2)
  56 {
  57   //---------------------------------------------------------------------
  58   // The main TPC tracker constructor
  59   //---------------------------------------------------------------------
  60   fInnerSec=new AliTPCSector[fkNIS];
  61   fOuterSec=new AliTPCSector[fkNOS];
  62
  63   Int_t i;
  64   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
  65   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
  66
  67   fN=0;  fSectors=0;
  68
  69   fClustersArray.Setup(par);
  70   fClustersArray.SetClusterType("AliTPCcluster");
  71   fClustersArray.ConnectTree("Segment Tree");
  72
  73   fSeeds=0;
  74 }
  75
  76 //_____________________________________________________________________________
  77 AliTPCtracker::~AliTPCtracker() {
  78   //------------------------------------------------------------------
  79   // TPC tracker destructor
  80   //------------------------------------------------------------------
  81   delete[] fInnerSec;
  82   delete[] fOuterSec;
  83   delete fSeeds;
  84 }
  85
  86 //_____________________________________________________________________________
  87 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
  88 {
  89   //
  90   // Parametrised error of the cluster reconstruction (pad direction)
  91   //
  92   // Sigma rphi
  93   const Float_t kArphi=0.41818e-2;
  94   const Float_t kBrphi=0.17460e-4;
  95   const Float_t kCrphi=0.30993e-2;
  96   const Float_t kDrphi=0.41061e-3;
  97
  98   pt=TMath::Abs(pt)*1000.;
  99   Double_t x=r/pt;
 100   tgl=TMath::Abs(tgl);
 101   Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
 102   if (s<0.4e-3) s=0.4e-3;
 103   s*=1.3; //Iouri Belikov
 104
 105   return s;
 106 }
 107
 108 //_____________________________________________________________________________
 109 Double_t SigmaZ2(Double_t r, Double_t tgl)
 110 {
 111   //
 112   // Parametrised error of the cluster reconstruction (drift direction)
 113   //
 114   // Sigma z
 115   const Float_t kAz=0.39614e-2;
 116   const Float_t kBz=0.22443e-4;
 117   const Float_t kCz=0.51504e-1;
 118
 119
 120   tgl=TMath::Abs(tgl);
 121   Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
 122   if (s<0.4e-3) s=0.4e-3;
 123   s*=1.3; //Iouri Belikov
 124
 125   return s;
 126 }
 127
 128 //_____________________________________________________________________________
 129 inline Double_t f1(Double_t x1,Double_t y1,
 130                    Double_t x2,Double_t y2,
 131                    Double_t x3,Double_t y3)
 132 {
 133   //-----------------------------------------------------------------
 134   // Initial approximation of the track curvature
 135   //-----------------------------------------------------------------
 136   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 137   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 138                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 139   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 140                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 141
 142   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 143
 144   return -xr*yr/sqrt(xr*xr+yr*yr);
 145 }
 146
 147
 148 //_____________________________________________________________________________
 149 inline Double_t f2(Double_t x1,Double_t y1,
 150                    Double_t x2,Double_t y2,
 151                    Double_t x3,Double_t y3)
 152 {
 153   //-----------------------------------------------------------------
 154   // Initial approximation of the track curvature times center of curvature
 155   //-----------------------------------------------------------------
 156   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 157   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 158                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 159   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 160                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 161
 162   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 163
 164   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 165 }
 166
 167 //_____________________________________________________________________________
 168 inline Double_t f3(Double_t x1,Double_t y1,
 169                    Double_t x2,Double_t y2,
 170                    Double_t z1,Double_t z2)
 171 {
 172   //-----------------------------------------------------------------
 173   // Initial approximation of the tangent of the track dip angle
 174   //-----------------------------------------------------------------
 175   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 176 }
 177
 178 //_____________________________________________________________________________
 179 void AliTPCtracker::LoadOuterSectors() {
 180   //-----------------------------------------------------------------
 181   // This function fills outer TPC sectors with clusters.
 182   //-----------------------------------------------------------------
 183   UInt_t index;
 184   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 185   for (Int_t i=0; i<j; i++) {
 186       AliSegmentID *s=fClustersArray.LoadEntry(i);
 187       Int_t sec,row;
 188       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 189       par->AdjustSectorRow(s->GetID(),sec,row);
 190       if (sec<fkNIS*2) continue;
 191       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 192       Int_t ncl=clrow->GetArray()->GetEntriesFast();
 193       while (ncl--) {
 194          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
 195          index=(((sec<<8)+row)<<16)+ncl;
 196          fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
 197       }
 198   }
 199
 200   fN=fkNOS;
 201   fSectors=fOuterSec;
 202 }
 203
 204 //_____________________________________________________________________________
 205 void AliTPCtracker::UnloadOuterSectors() {
 206   //-----------------------------------------------------------------
 207   // This function clears outer TPC sectors.
 208   //-----------------------------------------------------------------
 209   Int_t nup=fOuterSec->GetNRows();
 210   for (Int_t i=0; i<fkNOS; i++) {
 211     for (Int_t j=0; j<nup; j++) {
 212       if (fClustersArray.GetRow(i+fkNIS*2,j))
 213           fClustersArray.ClearRow(i+fkNIS*2,j);
 214       if (fClustersArray.GetRow(i+fkNIS*2+fkNOS,j))
 215           fClustersArray.ClearRow(i+fkNIS*2+fkNOS,j);
 216     }
 217   }
 218
 219   fN=0;
 220   fSectors=0;
 221 }
 222
 223 //_____________________________________________________________________________
 224 void AliTPCtracker::LoadInnerSectors() {
 225   //-----------------------------------------------------------------
 226   // This function fills inner TPC sectors with clusters.
 227   //-----------------------------------------------------------------
 228   UInt_t index;
 229   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 230   for (Int_t i=0; i<j; i++) {
 231       AliSegmentID *s=fClustersArray.LoadEntry(i);
 232       Int_t sec,row;
 233       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 234       par->AdjustSectorRow(s->GetID(),sec,row);
 235       if (sec>=fkNIS*2) continue;
 236       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 237       Int_t ncl=clrow->GetArray()->GetEntriesFast();
 238       while (ncl--) {
 239          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
 240          index=(((sec<<8)+row)<<16)+ncl;
 241          fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
 242       }
 243   }
 244
 245   fN=fkNIS;
 246   fSectors=fInnerSec;
 247 }
 248
 249 //_____________________________________________________________________________
 250 void AliTPCtracker::UnloadInnerSectors() {
 251   //-----------------------------------------------------------------
 252   // This function clears inner TPC sectors.
 253   //-----------------------------------------------------------------
 254   Int_t nlow=fInnerSec->GetNRows();
 255   for (Int_t i=0; i<fkNIS; i++) {
 256     for (Int_t j=0; j<nlow; j++) {
 257       if (fClustersArray.GetRow(i,j)) fClustersArray.ClearRow(i,j);
 258       if (fClustersArray.GetRow(i+fkNIS,j)) fClustersArray.ClearRow(i+fkNIS,j);
 259     }
 260   }
 261
 262   fN=0;
 263   fSectors=0;
 264 }
 265
 266 //_____________________________________________________________________________
 267 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
 268   //-----------------------------------------------------------------
 269   // This function tries to find a track prolongation.
 270   //-----------------------------------------------------------------
 271   Double_t xt=t.GetX();
 272   const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
 273                     Int_t(0.5*fSectors->GetNRows());
 274   Int_t tryAgain=kSKIP;
 275
 276   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
 277   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 278   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 279   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 280
 281   for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
 282     Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
 283     if (!t.PropagateTo(x)) return 0;
 284
 285     AliTPCcluster *cl=0;
 286     UInt_t index=0;
 287     Double_t maxchi2=kMaxCHI2;
 288     const AliTPCRow &krow=fSectors[s][nr];
 289     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
 290     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
 291     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
 292     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 293
 294     if (road>kMaxROAD) {
 295       if (t.GetNumberOfClusters()>4)
 296         cerr<<t.GetNumberOfClusters()
 297         <<"FindProlongation warning: Too broad road !\n";
 298       return 0;
 299     }
 300
 301     if (krow) {
 302       for (Int_t i=krow.Find(y-road); i<krow; i++) {
 303         AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 304         if (c->GetY() > y+road) break;
 305         if (c->IsUsed()) continue;
 306         if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
 307         Double_t chi2=t.GetPredictedChi2(c);
 308         if (chi2 > maxchi2) continue;
 309         maxchi2=chi2;
 310         cl=c;
 311         index=krow.GetIndex(i);
 312       }
 313     }
 314     if (cl) {
 315       Float_t l=fSectors->GetPadPitchWidth();
 316       t.SetSampledEdx(cl->GetQ()/l,t.GetNumberOfClusters());
 317       if (!t.Update(cl,maxchi2,index)) {
 318          if (!tryAgain--) return 0;
 319       } else tryAgain=kSKIP;
 320     } else {
 321       if (tryAgain==0) break;
 322       if (y > ymax) {
 323          s = (s+1) % fN;
 324          if (!t.Rotate(fSectors->GetAlpha())) return 0;
 325       } else if (y <-ymax) {
 326          s = (s-1+fN) % fN;
 327          if (!t.Rotate(-fSectors->GetAlpha())) return 0;
 328       }
 329       tryAgain--;
 330     }
 331   }
 332
 333   return 1;
 334 }
 335
 336 //_____________________________________________________________________________
 337 Int_t AliTPCtracker::FollowBackProlongation
 338 (AliTPCseed& seed, const AliTPCtrack &track) {
 339   //-----------------------------------------------------------------
 340   // This function propagates tracks back through the TPC
 341   //-----------------------------------------------------------------
 342   Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
 343   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 344   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 345   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 346
 347   Int_t idx=-1, sec=-1, row=-1;
 348   Int_t nc=seed.GetLabel(); //index of the cluster to start with
 349   if (nc--) {
 350      idx=track.GetClusterIndex(nc);
 351      sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 352   }
 353   if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
 354   else { if (sec <  2*fkNIS) row=-1; }
 355
 356   Int_t nr=fSectors->GetNRows();
 357   for (Int_t i=0; i<nr; i++) {
 358     Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
 359
 360     if (!seed.PropagateTo(x)) return 0;
 361
 362     Double_t y=seed.GetY();
 363     if (y > ymax) {
 364        s = (s+1) % fN;
 365        if (!seed.Rotate(fSectors->GetAlpha())) return 0;
 366     } else if (y <-ymax) {
 367        s = (s-1+fN) % fN;
 368        if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
 369     }
 370
 371     AliTPCcluster *cl=0;
 372     Int_t index=0;
 373     Double_t maxchi2=kMaxCHI2;
 374     Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
 375     Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
 376     Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
 377     Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
 378     if (road>kMaxROAD) {
 379       cerr<<seed.GetNumberOfClusters()
 380           <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
 381       return 0;
 382     }
 383
 384
 385     Int_t accepted=seed.GetNumberOfClusters();
 386     if (row==i) {
 387        //try to accept already found cluster
 388        AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
 389        Double_t chi2;
 390        if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
 391           index=idx; cl=c; maxchi2=chi2;
 392        } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
 393
 394        if (nc--) {
 395           idx=track.GetClusterIndex(nc);
 396           sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 397        }
 398        if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
 399        else { if (sec <  2*fkNIS) row=-1; }
 400
 401     }
 402     if (!cl) {
 403        //try to fill the gap
 404        const AliTPCRow &krow=fSectors[s][i];
 405        if (accepted>27)
 406        if (krow) {
 407           for (Int_t i=krow.Find(y-road); i<krow; i++) {
 408             AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 409             if (c->GetY() > y+road) break;
 410             if (c->IsUsed()) continue;
 411          if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
 412             Double_t chi2=seed.GetPredictedChi2(c);
 413             if (chi2 > maxchi2) continue;
 414             maxchi2=chi2;
 415             cl=c;
 416             index=krow.GetIndex(i);
 417           }
 418        }
 419     }
 420
 421     if (cl) {
 422       Float_t l=fSectors->GetPadPitchWidth();
 423       seed.SetSampledEdx(cl->GetQ()/l,seed.GetNumberOfClusters());
 424       seed.Update(cl,maxchi2,index);
 425     }
 426
 427   }
 428
 429   seed.SetLabel(nc);
 430
 431   return 1;
 432 }
 433
 434 //_____________________________________________________________________________
 435 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
 436   //-----------------------------------------------------------------
 437   // This function creates track seeds.
 438   //-----------------------------------------------------------------
 439   if (fSeeds==0) fSeeds=new TObjArray(15000);
 440
 441   Double_t x[5], c[15];
 442
 443   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
 444   Double_t cs=cos(alpha), sn=sin(alpha);
 445
 446   Double_t x1 =fOuterSec->GetX(i1);
 447   Double_t xx2=fOuterSec->GetX(i2);
 448
 449   for (Int_t ns=0; ns<fkNOS; ns++) {
 450     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
 451     Int_t nm=fOuterSec[ns][i2];
 452     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
 453     const AliTPCRow& kr1=fOuterSec[ns][i1];
 454     for (Int_t is=0; is < kr1; is++) {
 455       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
 456       for (Int_t js=0; js < nl+nm+nu; js++) {
 457         const AliTPCcluster *kcl;
 458         Double_t x2,   y2,   z2;
 459         Double_t x3=0.,y3=0.;
 460
 461         if (js<nl) {
 462           const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
 463           kcl=kr2[js];
 464           y2=kcl->GetY(); z2=kcl->GetZ();
 465           x2= xx2*cs+y2*sn;
 466           y2=-xx2*sn+y2*cs;
 467         } else
 468           if (js<nl+nm) {
 469             const AliTPCRow& kr2=fOuterSec[ns][i2];
 470             kcl=kr2[js-nl];
 471             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
 472           } else {
 473             const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
 474             kcl=kr2[js-nl-nm];
 475             y2=kcl->GetY(); z2=kcl->GetZ();
 476             x2=xx2*cs-y2*sn;
 477             y2=xx2*sn+y2*cs;
 478           }
 479
 480         Double_t zz=z1 - z1/x1*(x1-x2);
 481         if (TMath::Abs(zz-z2)>5.) continue;
 482
 483         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
 484         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
 485
 486         x[0]=y1;
 487         x[1]=z1;
 488         x[4]=f1(x1,y1,x2,y2,x3,y3);
 489         if (TMath::Abs(x[4]) >= 0.0066) continue;
 490         x[2]=f2(x1,y1,x2,y2,x3,y3);
 491         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
 492         x[3]=f3(x1,y1,x2,y2,z1,z2);
 493         if (TMath::Abs(x[3]) > 1.2) continue;
 494         Double_t a=asin(x[2]);
 495         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
 496         if (TMath::Abs(zv)>10.) continue;
 497
 498         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
 499         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
 500         Double_t sy3=100*0.025, sy=0.1, sz=0.1;
 501
 502         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
 503         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
 504         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
 505         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
 506         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
 507         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
 508         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
 509         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
 510         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
 511         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
 512
 513         c[0]=sy1;
 514         c[1]=0.;       c[2]=sz1;
 515         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
 516         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
 517                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
 518         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
 519         c[13]=f30*sy1*f40+f32*sy2*f42;
 520         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
 521
 522         UInt_t index=kr1.GetIndex(is);
 523         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
 524         Float_t l=fOuterSec->GetPadPitchWidth();
 525         track->SetSampledEdx(kr1[is]->GetQ()/l,0);
 526
 527         Int_t rc=FollowProlongation(*track, i2);
 528         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
 529         else fSeeds->AddLast(track);
 530       }
 531     }
 532   }
 533 }
 534
 535 //_____________________________________________________________________________
 536 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
 537   //-----------------------------------------------------------------
 538   // This function reades track seeds.
 539   //-----------------------------------------------------------------
 540   TDirectory *savedir=gDirectory;
 541
 542   TFile *in=(TFile*)inp;
 543   if (!in->IsOpen()) {
 544      cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
 545      return 1;
 546   }
 547
 548   in->cd();
 549   TTree *seedTree=(TTree*)in->Get("Seeds");
 550   if (!seedTree) {
 551      cerr<<"AliTPCtracker::ReadSeeds(): ";
 552      cerr<<"can't get a tree with track seeds !\n";
 553      return 2;
 554   }
 555   AliTPCtrack *seed=new AliTPCtrack;
 556   seedTree->SetBranchAddress("tracks",&seed);
 557
 558   if (fSeeds==0) fSeeds=new TObjArray(15000);
 559
 560   Int_t n=(Int_t)seedTree->GetEntries();
 561   for (Int_t i=0; i<n; i++) {
 562      seedTree->GetEvent(i);
 563      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
 564   }
 565
 566   delete seed;
 567   savedir->cd();
 568
 569   return 0;
 570 }
 571
 572 //_____________________________________________________________________________
 573 Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
 574   //-----------------------------------------------------------------
 575   // This is a track finder.
 576   //-----------------------------------------------------------------
 577   TDirectory *savedir=gDirectory;
 578
 579   if (inp) {
 580      TFile *in=(TFile*)inp;
 581      if (!in->IsOpen()) {
 582         cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
 583         return 1;
 584      }
 585   }
 586
 587   if (!out->IsOpen()) {
 588      cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
 589      return 2;
 590   }
 591
 592   out->cd();
 593   TTree tracktree("TPCf","Tree with TPC tracks");
 594   AliTPCtrack *iotrack=0;
 595   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 596
 597   LoadOuterSectors();
 598
 599   //find track seeds
 600   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 601   Int_t nrows=nlow+nup;
 602   if (fSeeds==0) {
 603      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
 604      MakeSeeds(nup-1, nup-1-gap);
 605      MakeSeeds(nup-1-shift, nup-1-shift-gap);
 606   }
 607   fSeeds->Sort();
 608
 609   //tracking in outer sectors
 610   Int_t nseed=fSeeds->GetEntriesFast();
 611   Int_t i;
 612   for (i=0; i<nseed; i++) {
 613     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 614     if (FollowProlongation(t)) {
 615        UseClusters(&t);
 616        continue;
 617     }
 618     delete fSeeds->RemoveAt(i);
 619   }
 620   UnloadOuterSectors();
 621
 622   //tracking in inner sectors
 623   LoadInnerSectors();
 624   Int_t found=0;
 625   for (i=0; i<nseed; i++) {
 626     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 627     if (!pt) continue;
 628     Int_t nc=t.GetNumberOfClusters();
 629
 630     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
 631     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 632     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 633     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
 634
 635     alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
 636
 637     if (t.Rotate(alpha)) {
 638        if (FollowProlongation(t)) {
 639           if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 640              t.CookdEdx();
 641              iotrack=pt;
 642              tracktree.Fill();
 643              UseClusters(&t,nc);
 644              cerr<<found++<<'\r';
 645           }
 646        }
 647     }
 648     delete fSeeds->RemoveAt(i);
 649   }
 650   UnloadInnerSectors();
 651
 652   tracktree.Write();
 653
 654   cerr<<"Number of found tracks : "<<found<<endl;
 655
 656   savedir->cd();
 657
 658   return 0;
 659 }
 660
 661 //_____________________________________________________________________________
 662 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
 663   //-----------------------------------------------------------------
 664   // This function propagates tracks back through the TPC.
 665   //-----------------------------------------------------------------
 666   fSeeds=new TObjArray(15000);
 667
 668   TFile *in=(TFile*)inp;
 669   TDirectory *savedir=gDirectory;
 670
 671   if (!in->IsOpen()) {
 672      cerr<<"AliTPCtracker::PropagateBack(): ";
 673      cerr<<"file with back propagated ITS tracks is not open !\n";
 674      return 1;
 675   }
 676
 677   if (!out->IsOpen()) {
 678      cerr<<"AliTPCtracker::PropagateBack(): ";
 679      cerr<<"file for back propagated TPC tracks is not open !\n";
 680      return 2;
 681   }
 682
 683   in->cd();
 684   TTree *bckTree=(TTree*)in->Get("ITSb");
 685   if (!bckTree) {
 686      cerr<<"AliTPCtracker::PropagateBack() ";
 687      cerr<<"can't get a tree with back propagated ITS tracks !\n";
 688      return 3;
 689   }
 690   AliTPCtrack *bckTrack=new AliTPCtrack;
 691   bckTree->SetBranchAddress("tracks",&bckTrack);
 692
 693   TTree *tpcTree=(TTree*)in->Get("TPCf");
 694   if (!tpcTree) {
 695      cerr<<"AliTPCtracker::PropagateBack() ";
 696      cerr<<"can't get a tree with TPC tracks !\n";
 697      return 4;
 698   }
 699   AliTPCtrack *tpcTrack=new AliTPCtrack;
 700   tpcTree->SetBranchAddress("tracks",&tpcTrack);
 701
 702 //*** Prepare an array of tracks to be back propagated
 703   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 704   Int_t nrows=nlow+nup;
 705
 706   TObjArray tracks(15000);
 707   Int_t i=0,j=0;
 708   Int_t tpcN=(Int_t)tpcTree->GetEntries();
 709   Int_t bckN=(Int_t)bckTree->GetEntries();
 710   if (j<bckN) bckTree->GetEvent(j++);
 711   for (i=0; i<tpcN; i++) {
 712      tpcTree->GetEvent(i);
 713      Double_t alpha=tpcTrack->GetAlpha();
 714      if (j<bckN &&
 715      TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
 716         if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
 717         fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
 718         bckTree->GetEvent(j++);
 719      } else {
 720         tpcTrack->ResetCovariance();
 721         fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
 722      }
 723      tracks.AddLast(new AliTPCtrack(*tpcTrack));
 724   }
 725
 726   out->cd();
 727   TTree backTree("TPCb","Tree with back propagated TPC tracks");
 728   AliTPCtrack *otrack=0;
 729   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 730
 731 //*** Back propagation through inner sectors
 732   LoadInnerSectors();
 733   Int_t nseed=fSeeds->GetEntriesFast();
 734   for (i=0; i<nseed; i++) {
 735      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
 736      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
 737
 738      Int_t nc=t.GetNumberOfClusters();
 739      s.SetLabel(nc-1); //set number of the cluster to start with
 740
 741      if (FollowBackProlongation(s,t)) {
 742         UseClusters(&s);
 743         continue;
 744      }
 745      delete fSeeds->RemoveAt(i);
 746   }
 747   UnloadInnerSectors();
 748
 749 //*** Back propagation through outer sectors
 750   LoadOuterSectors();
 751   Int_t found=0;
 752   for (i=0; i<nseed; i++) {
 753      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
 754      if (!ps) continue;
 755      Int_t nc=s.GetNumberOfClusters();
 756      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
 757
 758      Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
 759      if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 760      if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 761      Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
 762
 763      alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
 764      alpha-=s.GetAlpha();
 765
 766      if (s.Rotate(alpha)) {
 767         if (FollowBackProlongation(s,t)) {
 768            if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 769               s.CookdEdx();
 770               s.SetLabel(t.GetLabel());
 771               UseClusters(&s,nc);
 772               otrack=ps;
 773               backTree.Fill();
 774               cerr<<found++<<'\r';
 775               continue;
 776            }
 777         }
 778      }
 779      delete fSeeds->RemoveAt(i);
 780   }
 781   UnloadOuterSectors();
 782
 783   backTree.Write();
 784   savedir->cd();
 785   cerr<<"Number of seeds: "<<nseed<<endl;
 786   cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
 787   cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
 788
 789   delete bckTrack;
 790   delete tpcTrack;
 791
 792   return 0;
 793 }
 794
 795 //_________________________________________________________________________
 796 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
 797   //--------------------------------------------------------------------
 798   //       Return pointer to a given cluster
 799   //--------------------------------------------------------------------
 800   Int_t sec=(index&0xff000000)>>24;
 801   Int_t row=(index&0x00ff0000)>>16;
 802   Int_t ncl=(index&0x0000ffff)>>00;
 803
 804   AliTPCClustersRow *clrow=((AliTPCtracker*)this)->fClustersArray.GetRow(sec,row);;
 805   return (AliCluster*)(*clrow)[ncl];
 806 }
 807
 808 //__________________________________________________________________________
 809 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
 810   //--------------------------------------------------------------------
 811   //This function "cooks" a track label. If label<0, this track is fake.
 812   //--------------------------------------------------------------------
 813   Int_t noc=t->GetNumberOfClusters();
 814   Int_t *lb=new Int_t[noc];
 815   Int_t *mx=new Int_t[noc];
 816   AliCluster **clusters=new AliCluster*[noc];
 817
 818   Int_t i;
 819   for (i=0; i<noc; i++) {
 820      lb[i]=mx[i]=0;
 821      Int_t index=t->GetClusterIndex(i);
 822      clusters[i]=GetCluster(index);
 823   }
 824
 825   Int_t lab=123456789;
 826   for (i=0; i<noc; i++) {
 827     AliCluster *c=clusters[i];
 828     lab=TMath::Abs(c->GetLabel(0));
 829     Int_t j;
 830     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
 831     lb[j]=lab;
 832     (mx[j])++;
 833   }
 834
 835   Int_t max=0;
 836   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
 837
 838   for (i=0; i<noc; i++) {
 839     AliCluster *c=clusters[i];
 840     if (TMath::Abs(c->GetLabel(1)) == lab ||
 841         TMath::Abs(c->GetLabel(2)) == lab ) max++;
 842   }
 843
 844   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
 845
 846   else {
 847      Int_t tail=Int_t(0.10*noc);
 848      max=0;
 849      for (i=1; i<=tail; i++) {
 850        AliCluster *c=clusters[noc-i];
 851        if (lab == TMath::Abs(c->GetLabel(0)) ||
 852            lab == TMath::Abs(c->GetLabel(1)) ||
 853            lab == TMath::Abs(c->GetLabel(2))) max++;
 854      }
 855      if (max < Int_t(0.5*tail)) lab=-lab;
 856   }
 857
 858   t->SetLabel(lab);
 859
 860   delete[] lb;
 861   delete[] mx;
 862   delete[] clusters;
 863 }
 864
 865 //_________________________________________________________________________
 866 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
 867   //-----------------------------------------------------------------------
 868   // Setup inner sector
 869   //-----------------------------------------------------------------------
 870   if (f==0) {
 871      fAlpha=par->GetInnerAngle();
 872      fAlphaShift=par->GetInnerAngleShift();
 873      fPadPitchWidth=par->GetInnerPadPitchWidth();
 874      fPadPitchLength=par->GetInnerPadPitchLength();
 875
 876      fN=par->GetNRowLow();
 877      fRow=new AliTPCRow[fN];
 878      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
 879   } else {
 880      fAlpha=par->GetOuterAngle();
 881      fAlphaShift=par->GetOuterAngleShift();
 882      fPadPitchWidth=par->GetOuterPadPitchWidth();
 883      fPadPitchLength=par->GetOuterPadPitchLength();
 884
 885      fN=par->GetNRowUp();
 886      fRow=new AliTPCRow[fN];
 887      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiUp(i));
 888   }
 889 }
 890
 891 //_________________________________________________________________________
 892 void
 893 AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) {
 894   //-----------------------------------------------------------------------
 895   // Insert a cluster into this pad row in accordence with its y-coordinate
 896   //-----------------------------------------------------------------------
 897   if (fN==kMaxClusterPerRow) {
 898     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
 899   }
 900   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
 901   Int_t i=Find(c->GetY());
 902   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
 903   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
 904   fIndex[i]=index; fClusters[i]=c; fN++;
 905 }
 906
 907 //___________________________________________________________________
 908 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
 909   //-----------------------------------------------------------------------
 910   // Return the index of the nearest cluster
 911   //-----------------------------------------------------------------------
 912   if (y <= fClusters[0]->GetY()) return 0;
 913   if (y > fClusters[fN-1]->GetY()) return fN;
 914   Int_t b=0, e=fN-1, m=(b+e)/2;
 915   for (; b<e; m=(b+e)/2) {
 916     if (y > fClusters[m]->GetY()) b=m+1;
 917     else e=m;
 918   }
 919   return m;
 920 }
 921
 922 //_____________________________________________________________________________
 923 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
 924   //-----------------------------------------------------------------
 925   // This funtion calculates dE/dX within the "low" and "up" cuts.
 926   //-----------------------------------------------------------------
 927   Int_t i;
 928   Int_t nc=GetNumberOfClusters();
 929
 930   Int_t swap;//stupid sorting
 931   do {
 932     swap=0;
 933     for (i=0; i<nc-1; i++) {
 934       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
 935       Float_t tmp=fdEdxSample[i];
 936       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
 937       swap++;
 938     }
 939   } while (swap);
 940
 941   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
 942   Float_t dedx=0;
 943   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
 944   dedx /= (nu-nl+1);
 945   SetdEdx(dedx);
 946 }
 947
 948
 949