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