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