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