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