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