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.15  2001/11/08 16:36:33  hristov
  19 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)
  20
  21 Revision 1.14  2001/10/21 19:04:55  hristov
  22 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
  23
  24 Revision 1.13  2001/07/23 12:01:30  hristov
  25 Initialisation added
  26
  27 Revision 1.12  2001/07/20 14:32:44  kowal2
  28 Processing of many events possible now
  29
  30 Revision 1.11  2001/05/23 08:50:10  hristov
  31 Weird inline removed
  32
  33 Revision 1.10  2001/05/16 14:57:25  alibrary
  34 New files for folders and Stack
  35
  36 Revision 1.9  2001/05/11 07:16:56  hristov
  37 Fix needed on Sun and Alpha
  38
  39 Revision 1.8  2001/05/08 15:00:15  hristov
  40 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
  41
  42 Revision 1.5  2000/12/20 07:51:59  kowal2
  43 Changes suggested by Alessandra and Paolo to avoid overlapped
  44 data fields in encapsulated classes.
  45
  46 Revision 1.4  2000/11/02 07:27:16  kowal2
  47 code corrections
  48
  49 Revision 1.2  2000/06/30 12:07:50  kowal2
  50 Updated from the TPC-PreRelease branch
  51
  52 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
  53 Splitted from AliTPCtracking
  54
  55 */
  56
  57 //-------------------------------------------------------
  58 //          Implementation of the TPC tracker
  59 //
  60 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  61 //-------------------------------------------------------
  62
  63 #include <TObjArray.h>
  64 #include <TFile.h>
  65 #include <TTree.h>
  66 #include <iostream.h>
  67
  68 #include "AliTPCtracker.h"
  69 #include "AliTPCcluster.h"
  70 #include "AliTPCParam.h"
  71 #include "AliTPCClustersRow.h"
  72
  73 //_____________________________________________________________________________
  74 AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn):
  75 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
  76 {
  77   //---------------------------------------------------------------------
  78   // The main TPC tracker constructor
  79   //---------------------------------------------------------------------
  80   fInnerSec=new AliTPCSector[fkNIS];
  81   fOuterSec=new AliTPCSector[fkNOS];
  82
  83   Int_t i;
  84   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
  85   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
  86
  87   fN=0;  fSectors=0;
  88
  89   fClustersArray.Setup(par);
  90   fClustersArray.SetClusterType("AliTPCcluster");
  91
  92   char   cname[100];
  93   if (eventn==-1) {
  94     sprintf(cname,"TreeC_TPC");
  95   }
  96   else {
  97     sprintf(cname,"TreeC_TPC_%d",eventn);
  98   }
  99
 100   fClustersArray.ConnectTree(cname);
 101
 102   fEventN = eventn;
 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::LoadOuterSectors() {
 213   //-----------------------------------------------------------------
 214   // This function fills outer TPC sectors with clusters.
 215   //-----------------------------------------------------------------
 216   UInt_t index;
 217   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 218   for (Int_t i=0; i<j; i++) {
 219       AliSegmentID *s=fClustersArray.LoadEntry(i);
 220       Int_t sec,row;
 221       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 222       par->AdjustSectorRow(s->GetID(),sec,row);
 223       if (sec<fkNIS*2) continue;
 224       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 225       Int_t ncl=clrow->GetArray()->GetEntriesFast();
 226       while (ncl--) {
 227          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
 228          index=(((sec<<8)+row)<<16)+ncl;
 229          fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
 230       }
 231   }
 232
 233   fN=fkNOS;
 234   fSectors=fOuterSec;
 235 }
 236
 237 //_____________________________________________________________________________
 238 void AliTPCtracker::UnloadOuterSectors() {
 239   //-----------------------------------------------------------------
 240   // This function clears outer TPC sectors.
 241   //-----------------------------------------------------------------
 242   Int_t nup=fOuterSec->GetNRows();
 243   for (Int_t i=0; i<fkNOS; i++) {
 244     for (Int_t j=0; j<nup; j++) {
 245       if (fClustersArray.GetRow(i+fkNIS*2,j))
 246           fClustersArray.ClearRow(i+fkNIS*2,j);
 247       if (fClustersArray.GetRow(i+fkNIS*2+fkNOS,j))
 248           fClustersArray.ClearRow(i+fkNIS*2+fkNOS,j);
 249     }
 250   }
 251
 252   fN=0;
 253   fSectors=0;
 254 }
 255
 256 //_____________________________________________________________________________
 257 void AliTPCtracker::LoadInnerSectors() {
 258   //-----------------------------------------------------------------
 259   // This function fills inner TPC sectors with clusters.
 260   //-----------------------------------------------------------------
 261   UInt_t index;
 262   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 263   for (Int_t i=0; i<j; i++) {
 264       AliSegmentID *s=fClustersArray.LoadEntry(i);
 265       Int_t sec,row;
 266       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 267       par->AdjustSectorRow(s->GetID(),sec,row);
 268       if (sec>=fkNIS*2) continue;
 269       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 270       Int_t ncl=clrow->GetArray()->GetEntriesFast();
 271       while (ncl--) {
 272          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
 273          index=(((sec<<8)+row)<<16)+ncl;
 274          fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
 275       }
 276   }
 277
 278   fN=fkNIS;
 279   fSectors=fInnerSec;
 280 }
 281
 282 //_____________________________________________________________________________
 283 void AliTPCtracker::UnloadInnerSectors() {
 284   //-----------------------------------------------------------------
 285   // This function clears inner TPC sectors.
 286   //-----------------------------------------------------------------
 287   Int_t nlow=fInnerSec->GetNRows();
 288   for (Int_t i=0; i<fkNIS; i++) {
 289     for (Int_t j=0; j<nlow; j++) {
 290       if (fClustersArray.GetRow(i,j)) fClustersArray.ClearRow(i,j);
 291       if (fClustersArray.GetRow(i+fkNIS,j)) fClustersArray.ClearRow(i+fkNIS,j);
 292     }
 293   }
 294
 295   fN=0;
 296   fSectors=0;
 297 }
 298
 299 //_____________________________________________________________________________
 300 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
 301   //-----------------------------------------------------------------
 302   // This function tries to find a track prolongation.
 303   //-----------------------------------------------------------------
 304   Double_t xt=t.GetX();
 305   const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
 306                     Int_t(0.5*fSectors->GetNRows());
 307   Int_t tryAgain=kSKIP;
 308
 309   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
 310   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 311   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 312   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 313
 314   for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
 315     Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
 316     if (!t.PropagateTo(x)) return 0;
 317
 318     AliTPCcluster *cl=0;
 319     UInt_t index=0;
 320     Double_t maxchi2=kMaxCHI2;
 321     const AliTPCRow &krow=fSectors[s][nr];
 322     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
 323     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
 324     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
 325     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 326
 327     if (road>kMaxROAD) {
 328       if (t.GetNumberOfClusters()>4)
 329         cerr<<t.GetNumberOfClusters()
 330         <<"FindProlongation warning: Too broad road !\n";
 331       return 0;
 332     }
 333
 334     if (krow) {
 335       for (Int_t i=krow.Find(y-road); i<krow; i++) {
 336         AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 337         if (c->GetY() > y+road) break;
 338         if (c->IsUsed()) continue;
 339         if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
 340         Double_t chi2=t.GetPredictedChi2(c);
 341         if (chi2 > maxchi2) continue;
 342         maxchi2=chi2;
 343         cl=c;
 344         index=krow.GetIndex(i);
 345       }
 346     }
 347     if (cl) {
 348       Float_t l=fSectors->GetPadPitchWidth();
 349       t.SetSampledEdx(cl->GetQ()/l,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 >= 2*fkNIS) row=-1; }
 387   else { if (sec <  2*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 >= 2*fkNIS) row=-1; }
 432        else { if (sec <  2*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       seed.SetSampledEdx(cl->GetQ()/l,seed.GetNumberOfClusters());
 457       seed.Update(cl,maxchi2,index);
 458     }
 459
 460   }
 461
 462   seed.SetLabel(nc);
 463
 464   return 1;
 465 }
 466
 467 //_____________________________________________________________________________
 468 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
 469   //-----------------------------------------------------------------
 470   // This function creates track seeds.
 471   //-----------------------------------------------------------------
 472   if (fSeeds==0) fSeeds=new TObjArray(15000);
 473
 474   Double_t x[5], c[15];
 475
 476   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
 477   Double_t cs=cos(alpha), sn=sin(alpha);
 478
 479   Double_t x1 =fOuterSec->GetX(i1);
 480   Double_t xx2=fOuterSec->GetX(i2);
 481
 482   for (Int_t ns=0; ns<fkNOS; ns++) {
 483     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
 484     Int_t nm=fOuterSec[ns][i2];
 485     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
 486     const AliTPCRow& kr1=fOuterSec[ns][i1];
 487     for (Int_t is=0; is < kr1; is++) {
 488       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
 489       for (Int_t js=0; js < nl+nm+nu; js++) {
 490         const AliTPCcluster *kcl;
 491         Double_t x2,   y2,   z2;
 492         Double_t x3=GetX(), y3=GetY(), z3=GetZ();
 493
 494         if (js<nl) {
 495           const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
 496           kcl=kr2[js];
 497           y2=kcl->GetY(); z2=kcl->GetZ();
 498           x2= xx2*cs+y2*sn;
 499           y2=-xx2*sn+y2*cs;
 500         } else
 501           if (js<nl+nm) {
 502             const AliTPCRow& kr2=fOuterSec[ns][i2];
 503             kcl=kr2[js-nl];
 504             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
 505           } else {
 506             const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
 507             kcl=kr2[js-nl-nm];
 508             y2=kcl->GetY(); z2=kcl->GetZ();
 509             x2=xx2*cs-y2*sn;
 510             y2=xx2*sn+y2*cs;
 511           }
 512
 513         Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
 514         if (TMath::Abs(zz-z2)>5.) continue;
 515
 516         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
 517         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
 518
 519         x[0]=y1;
 520         x[1]=z1;
 521         x[4]=f1(x1,y1,x2,y2,x3,y3);
 522         if (TMath::Abs(x[4]) >= 0.0066) continue;
 523         x[2]=f2(x1,y1,x2,y2,x3,y3);
 524         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
 525         x[3]=f3(x1,y1,x2,y2,z1,z2);
 526         if (TMath::Abs(x[3]) > 1.2) continue;
 527         Double_t a=asin(x[2]);
 528         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
 529         if (TMath::Abs(zv-z3)>10.) continue;
 530
 531         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
 532         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
 533         Double_t sy3=400*3./12., 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=fOuterSec->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   out->cd();
 629
 630   char   tname[100];
 631   sprintf(tname,"TreeT_TPC_%d",fEventN);
 632   TTree tracktree(tname,"Tree with TPC tracks");
 633   AliTPCtrack *iotrack=0;
 634   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 635
 636   LoadOuterSectors();
 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      MakeSeeds(nup-1, nup-1-gap);
 644      MakeSeeds(nup-1-shift, nup-1-shift-gap);
 645   }
 646   fSeeds->Sort();
 647
 648   //tracking in outer sectors
 649   Int_t nseed=fSeeds->GetEntriesFast();
 650   Int_t i;
 651   for (i=0; i<nseed; i++) {
 652     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 653     if (FollowProlongation(t)) {
 654        UseClusters(&t);
 655        continue;
 656     }
 657     delete fSeeds->RemoveAt(i);
 658   }
 659   //UnloadOuterSectors();
 660
 661   //tracking in inner sectors
 662   LoadInnerSectors();
 663   Int_t found=0;
 664   for (i=0; i<nseed; i++) {
 665     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 666     if (!pt) continue;
 667     Int_t nc=t.GetNumberOfClusters();
 668
 669     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
 670     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 671     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 672     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
 673
 674     alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
 675
 676     if (t.Rotate(alpha)) {
 677        if (FollowProlongation(t)) {
 678           if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 679              t.CookdEdx();
 680              CookLabel(pt,0.1); //For comparison only
 681              iotrack=pt;
 682              tracktree.Fill();
 683              UseClusters(&t,nc);
 684              cerr<<found++<<'\r';
 685           }
 686        }
 687     }
 688     delete fSeeds->RemoveAt(i);
 689   }
 690   UnloadInnerSectors();
 691   UnloadOuterSectors();
 692
 693   tracktree.Write();
 694
 695   cerr<<"Number of found tracks : "<<found<<endl;
 696
 697   savedir->cd();
 698
 699   return 0;
 700 }
 701
 702 //_____________________________________________________________________________
 703 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
 704   //-----------------------------------------------------------------
 705   // This function propagates tracks back through the TPC.
 706   //-----------------------------------------------------------------
 707   fSeeds=new TObjArray(15000);
 708
 709   TFile *in=(TFile*)inp;
 710   TDirectory *savedir=gDirectory;
 711
 712   if (!in->IsOpen()) {
 713      cerr<<"AliTPCtracker::PropagateBack(): ";
 714      cerr<<"file with back propagated ITS tracks is not open !\n";
 715      return 1;
 716   }
 717
 718   if (!out->IsOpen()) {
 719      cerr<<"AliTPCtracker::PropagateBack(): ";
 720      cerr<<"file for back propagated TPC tracks is not open !\n";
 721      return 2;
 722   }
 723
 724   in->cd();
 725   TTree *bckTree=(TTree*)in->Get("TreeT_ITSb_0");
 726   if (!bckTree) {
 727      cerr<<"AliTPCtracker::PropagateBack() ";
 728      cerr<<"can't get a tree with back propagated ITS tracks !\n";
 729      return 3;
 730   }
 731   AliTPCtrack *bckTrack=new AliTPCtrack;
 732   bckTree->SetBranchAddress("tracks",&bckTrack);
 733
 734   TTree *tpcTree=(TTree*)in->Get("TreeT_TPC_0");
 735   if (!tpcTree) {
 736      cerr<<"AliTPCtracker::PropagateBack() ";
 737      cerr<<"can't get a tree with TPC tracks !\n";
 738      return 4;
 739   }
 740   AliTPCtrack *tpcTrack=new AliTPCtrack;
 741   tpcTree->SetBranchAddress("tracks",&tpcTrack);
 742
 743 //*** Prepare an array of tracks to be back propagated
 744   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 745   Int_t nrows=nlow+nup;
 746
 747   TObjArray tracks(15000);
 748   Int_t i=0,j=0;
 749   Int_t tpcN=(Int_t)tpcTree->GetEntries();
 750   Int_t bckN=(Int_t)bckTree->GetEntries();
 751   if (j<bckN) bckTree->GetEvent(j++);
 752   for (i=0; i<tpcN; i++) {
 753      tpcTree->GetEvent(i);
 754      Double_t alpha=tpcTrack->GetAlpha();
 755      if (j<bckN &&
 756      TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
 757         if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
 758         fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
 759         bckTree->GetEvent(j++);
 760      } else {
 761         tpcTrack->ResetCovariance();
 762         fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
 763      }
 764      tracks.AddLast(new AliTPCtrack(*tpcTrack));
 765   }
 766
 767   out->cd();
 768   TTree backTree("TreeT_TPCb_0","Tree with back propagated TPC tracks");
 769   AliTPCtrack *otrack=0;
 770   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 771
 772 //*** Back propagation through inner sectors
 773   LoadInnerSectors();
 774   Int_t nseed=fSeeds->GetEntriesFast();
 775   for (i=0; i<nseed; i++) {
 776      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
 777      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
 778
 779      Int_t nc=t.GetNumberOfClusters();
 780      s.SetLabel(nc-1); //set number of the cluster to start with
 781
 782      if (FollowBackProlongation(s,t)) {
 783         UseClusters(&s);
 784         continue;
 785      }
 786      delete fSeeds->RemoveAt(i);
 787   }
 788   UnloadInnerSectors();
 789
 790 //*** Back propagation through outer sectors
 791   LoadOuterSectors();
 792   Int_t found=0;
 793   for (i=0; i<nseed; i++) {
 794      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
 795      if (!ps) continue;
 796      Int_t nc=s.GetNumberOfClusters();
 797      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
 798
 799      Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
 800      if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 801      if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 802      Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
 803
 804      alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
 805      alpha-=s.GetAlpha();
 806
 807      if (s.Rotate(alpha)) {
 808         if (FollowBackProlongation(s,t)) {
 809            if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 810               s.CookdEdx();
 811               s.SetLabel(t.GetLabel());
 812               UseClusters(&s,nc);
 813               otrack=ps;
 814               backTree.Fill();
 815               cerr<<found++<<'\r';
 816               continue;
 817            }
 818         }
 819      }
 820      delete fSeeds->RemoveAt(i);
 821   }
 822   UnloadOuterSectors();
 823
 824   backTree.Write();
 825   savedir->cd();
 826   cerr<<"Number of seeds: "<<nseed<<endl;
 827   cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
 828   cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
 829
 830   delete bckTrack;
 831   delete tpcTrack;
 832
 833   delete bckTree; //Thanks to Mariana Bondila
 834   delete tpcTree; //Thanks to Mariana Bondila
 835
 836   return 0;
 837 }
 838
 839 //_________________________________________________________________________
 840 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
 841   //--------------------------------------------------------------------
 842   //       Return pointer to a given cluster
 843   //--------------------------------------------------------------------
 844   Int_t sec=(index&0xff000000)>>24;
 845   Int_t row=(index&0x00ff0000)>>16;
 846   Int_t ncl=(index&0x0000ffff)>>00;
 847
 848   AliTPCClustersRow *clrow=((AliTPCtracker *) this)->fClustersArray.GetRow(sec,row);
 849   return (AliCluster*)(*clrow)[ncl];
 850 }
 851
 852 //__________________________________________________________________________
 853 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
 854   //--------------------------------------------------------------------
 855   //This function "cooks" a track label. If label<0, this track is fake.
 856   //--------------------------------------------------------------------
 857   Int_t noc=t->GetNumberOfClusters();
 858   Int_t *lb=new Int_t[noc];
 859   Int_t *mx=new Int_t[noc];
 860   AliCluster **clusters=new AliCluster*[noc];
 861
 862   Int_t i;
 863   for (i=0; i<noc; i++) {
 864      lb[i]=mx[i]=0;
 865      Int_t index=t->GetClusterIndex(i);
 866      clusters[i]=GetCluster(index);
 867   }
 868
 869   Int_t lab=123456789;
 870   for (i=0; i<noc; i++) {
 871     AliCluster *c=clusters[i];
 872     lab=TMath::Abs(c->GetLabel(0));
 873     Int_t j;
 874     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
 875     lb[j]=lab;
 876     (mx[j])++;
 877   }
 878
 879   Int_t max=0;
 880   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
 881
 882   for (i=0; i<noc; i++) {
 883     AliCluster *c=clusters[i];
 884     if (TMath::Abs(c->GetLabel(1)) == lab ||
 885         TMath::Abs(c->GetLabel(2)) == lab ) max++;
 886   }
 887
 888   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
 889
 890   else {
 891      Int_t tail=Int_t(0.10*noc);
 892      max=0;
 893      for (i=1; i<=tail; i++) {
 894        AliCluster *c=clusters[noc-i];
 895        if (lab == TMath::Abs(c->GetLabel(0)) ||
 896            lab == TMath::Abs(c->GetLabel(1)) ||
 897            lab == TMath::Abs(c->GetLabel(2))) max++;
 898      }
 899      if (max < Int_t(0.5*tail)) lab=-lab;
 900   }
 901
 902   t->SetLabel(lab);
 903
 904   delete[] lb;
 905   delete[] mx;
 906   delete[] clusters;
 907 }
 908
 909 //_________________________________________________________________________
 910 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
 911   //-----------------------------------------------------------------------
 912   // Setup inner sector
 913   //-----------------------------------------------------------------------
 914   if (f==0) {
 915      fAlpha=par->GetInnerAngle();
 916      fAlphaShift=par->GetInnerAngleShift();
 917      fPadPitchWidth=par->GetInnerPadPitchWidth();
 918      fPadPitchLength=par->GetInnerPadPitchLength();
 919
 920      fN=par->GetNRowLow();
 921      fRow=new AliTPCRow[fN];
 922      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
 923   } else {
 924      fAlpha=par->GetOuterAngle();
 925      fAlphaShift=par->GetOuterAngleShift();
 926      fPadPitchWidth=par->GetOuterPadPitchWidth();
 927      fPadPitchLength=par->GetOuterPadPitchLength();
 928
 929      fN=par->GetNRowUp();
 930      fRow=new AliTPCRow[fN];
 931      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiUp(i));
 932   }
 933 }
 934
 935 //_________________________________________________________________________
 936 void
 937 AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) {
 938   //-----------------------------------------------------------------------
 939   // Insert a cluster into this pad row in accordence with its y-coordinate
 940   //-----------------------------------------------------------------------
 941   if (fN==kMaxClusterPerRow) {
 942     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
 943   }
 944   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
 945   Int_t i=Find(c->GetY());
 946   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
 947   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
 948   fIndex[i]=index; fClusters[i]=c; fN++;
 949 }
 950
 951 //___________________________________________________________________
 952 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
 953   //-----------------------------------------------------------------------
 954   // Return the index of the nearest cluster
 955   //-----------------------------------------------------------------------
 956   if (y <= fClusters[0]->GetY()) return 0;
 957   if (y > fClusters[fN-1]->GetY()) return fN;
 958   Int_t b=0, e=fN-1, m=(b+e)/2;
 959   for (; b<e; m=(b+e)/2) {
 960     if (y > fClusters[m]->GetY()) b=m+1;
 961     else e=m;
 962   }
 963   return m;
 964 }
 965
 966 //_____________________________________________________________________________
 967 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
 968   //-----------------------------------------------------------------
 969   // This funtion calculates dE/dX within the "low" and "up" cuts.
 970   //-----------------------------------------------------------------
 971   Int_t i;
 972   Int_t nc=GetNumberOfClusters();
 973
 974   Int_t swap;//stupid sorting
 975   do {
 976     swap=0;
 977     for (i=0; i<nc-1; i++) {
 978       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
 979       Float_t tmp=fdEdxSample[i];
 980       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
 981       swap++;
 982     }
 983   } while (swap);
 984
 985   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
 986   Float_t dedx=0;
 987   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
 988   dedx /= (nu-nl+1);
 989   SetdEdx(dedx);
 990
 991   //Very rough PID
 992   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
 993
 994   if (p<0.6) {
 995     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
 996     if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
 997     SetMass(0.93827); return;
 998   }
 999
1000   if (p<1.2) {
1001     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1002     SetMass(0.93827); return;
1003   }
1004
1005   SetMass(0.13957); return;
1006
1007 }
1008
1009
1010