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