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