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