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