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