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