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