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.16  2001/11/08 16:39:03  hristov
  19 Additional protection (M.Masera)
  20
  21 Revision 1.15  2001/11/08 16:36:33  hristov
  22 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)
  23
  24 Revision 1.14  2001/10/21 19:04:55  hristov
  25 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
  26
  27 Revision 1.13  2001/07/23 12:01:30  hristov
  28 Initialisation added
  29
  30 Revision 1.12  2001/07/20 14:32:44  kowal2
  31 Processing of many events possible now
  32
  33 Revision 1.11  2001/05/23 08:50:10  hristov
  34 Weird inline removed
  35
  36 Revision 1.10  2001/05/16 14:57:25  alibrary
  37 New files for folders and Stack
  38
  39 Revision 1.9  2001/05/11 07:16:56  hristov
  40 Fix needed on Sun and Alpha
  41
  42 Revision 1.8  2001/05/08 15:00:15  hristov
  43 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
  44
  45 Revision 1.5  2000/12/20 07:51:59  kowal2
  46 Changes suggested by Alessandra and Paolo to avoid overlapped
  47 data fields in encapsulated classes.
  48
  49 Revision 1.4  2000/11/02 07:27:16  kowal2
  50 code corrections
  51
  52 Revision 1.2  2000/06/30 12:07:50  kowal2
  53 Updated from the TPC-PreRelease branch
  54
  55 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
  56 Splitted from AliTPCtracking
  57
  58 */
  59
  60 //-------------------------------------------------------
  61 //          Implementation of the TPC tracker
  62 //
  63 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  64 //-------------------------------------------------------
  65
  66 #include <TObjArray.h>
  67 #include <TFile.h>
  68 #include <TTree.h>
  69 #include <iostream.h>
  70
  71 #include "AliTPCtracker.h"
  72 #include "AliTPCcluster.h"
  73 #include "AliTPCParam.h"
  74 #include "AliTPCClustersRow.h"
  75
  76 //_____________________________________________________________________________
  77 AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn):
  78 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
  79 {
  80   //---------------------------------------------------------------------
  81   // The main TPC tracker constructor
  82   //---------------------------------------------------------------------
  83   fInnerSec=new AliTPCSector[fkNIS];
  84   fOuterSec=new AliTPCSector[fkNOS];
  85
  86   Int_t i;
  87   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
  88   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
  89
  90   fN=0;  fSectors=0;
  91
  92   fClustersArray.Setup(par);
  93   fClustersArray.SetClusterType("AliTPCcluster");
  94
  95   char   cname[100];
  96   if (eventn==-1) {
  97     sprintf(cname,"TreeC_TPC");
  98   }
  99   else {
 100     sprintf(cname,"TreeC_TPC_%d",eventn);
 101   }
 102
 103   fClustersArray.ConnectTree(cname);
 104
 105   fEventN = eventn;
 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::LoadOuterSectors() {
 216   //-----------------------------------------------------------------
 217   // This function fills outer TPC sectors with clusters.
 218   //-----------------------------------------------------------------
 219   UInt_t index;
 220   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 221   for (Int_t i=0; i<j; i++) {
 222       AliSegmentID *s=fClustersArray.LoadEntry(i);
 223       Int_t sec,row;
 224       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 225       par->AdjustSectorRow(s->GetID(),sec,row);
 226       if (sec<fkNIS*2) continue;
 227       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 228       Int_t ncl=clrow->GetArray()->GetEntriesFast();
 229       while (ncl--) {
 230          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
 231          index=(((sec<<8)+row)<<16)+ncl;
 232          fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
 233       }
 234   }
 235
 236   fN=fkNOS;
 237   fSectors=fOuterSec;
 238 }
 239
 240 //_____________________________________________________________________________
 241 void AliTPCtracker::UnloadOuterSectors() {
 242   //-----------------------------------------------------------------
 243   // This function clears outer TPC sectors.
 244   //-----------------------------------------------------------------
 245   Int_t nup=fOuterSec->GetNRows();
 246   for (Int_t i=0; i<fkNOS; i++) {
 247     for (Int_t j=0; j<nup; j++) {
 248       if (fClustersArray.GetRow(i+fkNIS*2,j))
 249           fClustersArray.ClearRow(i+fkNIS*2,j);
 250       if (fClustersArray.GetRow(i+fkNIS*2+fkNOS,j))
 251           fClustersArray.ClearRow(i+fkNIS*2+fkNOS,j);
 252     }
 253   }
 254
 255   fN=0;
 256   fSectors=0;
 257 }
 258
 259 //_____________________________________________________________________________
 260 void AliTPCtracker::LoadInnerSectors() {
 261   //-----------------------------------------------------------------
 262   // This function fills inner TPC sectors with clusters.
 263   //-----------------------------------------------------------------
 264   UInt_t index;
 265   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 266   for (Int_t i=0; i<j; i++) {
 267       AliSegmentID *s=fClustersArray.LoadEntry(i);
 268       Int_t sec,row;
 269       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 270       par->AdjustSectorRow(s->GetID(),sec,row);
 271       if (sec>=fkNIS*2) continue;
 272       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 273       Int_t ncl=clrow->GetArray()->GetEntriesFast();
 274       while (ncl--) {
 275          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
 276          index=(((sec<<8)+row)<<16)+ncl;
 277          fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
 278       }
 279   }
 280
 281   fN=fkNIS;
 282   fSectors=fInnerSec;
 283 }
 284
 285 //_____________________________________________________________________________
 286 void AliTPCtracker::UnloadInnerSectors() {
 287   //-----------------------------------------------------------------
 288   // This function clears inner TPC sectors.
 289   //-----------------------------------------------------------------
 290   Int_t nlow=fInnerSec->GetNRows();
 291   for (Int_t i=0; i<fkNIS; i++) {
 292     for (Int_t j=0; j<nlow; j++) {
 293       if (fClustersArray.GetRow(i,j)) fClustersArray.ClearRow(i,j);
 294       if (fClustersArray.GetRow(i+fkNIS,j)) fClustersArray.ClearRow(i+fkNIS,j);
 295     }
 296   }
 297
 298   fN=0;
 299   fSectors=0;
 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       t.SetSampledEdx(cl->GetQ()/l,t.GetNumberOfClusters());
 354       if (!t.Update(cl,maxchi2,index)) {
 355          if (!tryAgain--) return 0;
 356       } else tryAgain=kSKIP;
 357     } else {
 358       if (tryAgain==0) break;
 359       if (y > ymax) {
 360          s = (s+1) % fN;
 361          if (!t.Rotate(fSectors->GetAlpha())) return 0;
 362       } else if (y <-ymax) {
 363          s = (s-1+fN) % fN;
 364          if (!t.Rotate(-fSectors->GetAlpha())) return 0;
 365       }
 366       tryAgain--;
 367     }
 368   }
 369
 370   return 1;
 371 }
 372
 373 //_____________________________________________________________________________
 374 Int_t AliTPCtracker::FollowBackProlongation
 375 (AliTPCseed& seed, const AliTPCtrack &track) {
 376   //-----------------------------------------------------------------
 377   // This function propagates tracks back through the TPC
 378   //-----------------------------------------------------------------
 379   Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
 380   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 381   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 382   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 383
 384   Int_t idx=-1, sec=-1, row=-1;
 385   Int_t nc=seed.GetLabel(); //index of the cluster to start with
 386   if (nc--) {
 387      idx=track.GetClusterIndex(nc);
 388      sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 389   }
 390   if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
 391   else { if (sec <  2*fkNIS) row=-1; }
 392
 393   Int_t nr=fSectors->GetNRows();
 394   for (Int_t i=0; i<nr; i++) {
 395     Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
 396
 397     if (!seed.PropagateTo(x)) return 0;
 398
 399     Double_t y=seed.GetY();
 400     if (y > ymax) {
 401        s = (s+1) % fN;
 402        if (!seed.Rotate(fSectors->GetAlpha())) return 0;
 403     } else if (y <-ymax) {
 404        s = (s-1+fN) % fN;
 405        if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
 406     }
 407
 408     AliTPCcluster *cl=0;
 409     Int_t index=0;
 410     Double_t maxchi2=kMaxCHI2;
 411     Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
 412     Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
 413     Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
 414     Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
 415     if (road>kMaxROAD) {
 416       cerr<<seed.GetNumberOfClusters()
 417           <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
 418       return 0;
 419     }
 420
 421
 422     Int_t accepted=seed.GetNumberOfClusters();
 423     if (row==i) {
 424        //try to accept already found cluster
 425        AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
 426        Double_t chi2;
 427        if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
 428           index=idx; cl=c; maxchi2=chi2;
 429        } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
 430
 431        if (nc--) {
 432           idx=track.GetClusterIndex(nc);
 433           sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 434        }
 435        if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
 436        else { if (sec <  2*fkNIS) row=-1; }
 437
 438     }
 439     if (!cl) {
 440        //try to fill the gap
 441        const AliTPCRow &krow=fSectors[s][i];
 442        if (accepted>27)
 443        if (krow) {
 444           for (Int_t i=krow.Find(y-road); i<krow; i++) {
 445             AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 446             if (c->GetY() > y+road) break;
 447             if (c->IsUsed()) continue;
 448          if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
 449             Double_t chi2=seed.GetPredictedChi2(c);
 450             if (chi2 > maxchi2) continue;
 451             maxchi2=chi2;
 452             cl=c;
 453             index=krow.GetIndex(i);
 454           }
 455        }
 456     }
 457
 458     if (cl) {
 459       Float_t l=fSectors->GetPadPitchWidth();
 460       seed.SetSampledEdx(cl->GetQ()/l,seed.GetNumberOfClusters());
 461       seed.Update(cl,maxchi2,index);
 462     }
 463
 464   }
 465
 466   seed.SetLabel(nc);
 467
 468   return 1;
 469 }
 470
 471 //_____________________________________________________________________________
 472 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
 473   //-----------------------------------------------------------------
 474   // This function creates track seeds.
 475   //-----------------------------------------------------------------
 476   if (fSeeds==0) fSeeds=new TObjArray(15000);
 477
 478   Double_t x[5], c[15];
 479
 480   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
 481   Double_t cs=cos(alpha), sn=sin(alpha);
 482
 483   Double_t x1 =fOuterSec->GetX(i1);
 484   Double_t xx2=fOuterSec->GetX(i2);
 485
 486   for (Int_t ns=0; ns<fkNOS; ns++) {
 487     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
 488     Int_t nm=fOuterSec[ns][i2];
 489     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
 490     const AliTPCRow& kr1=fOuterSec[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=fOuterSec[(ns-1+fkNOS)%fkNOS][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=fOuterSec[ns][i2];
 507             kcl=kr2[js-nl];
 508             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
 509           } else {
 510             const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][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
 539         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
 540         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
 541         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
 542         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
 543         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
 544         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
 545         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
 546         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
 547         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
 548         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
 549
 550         c[0]=sy1;
 551         c[1]=0.;       c[2]=sz1;
 552         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
 553         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
 554                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
 555         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
 556         c[13]=f30*sy1*f40+f32*sy2*f42;
 557         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
 558
 559         UInt_t index=kr1.GetIndex(is);
 560         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
 561         Float_t l=fOuterSec->GetPadPitchWidth();
 562         track->SetSampledEdx(kr1[is]->GetQ()/l,0);
 563
 564         Int_t rc=FollowProlongation(*track, i2);
 565         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
 566         else fSeeds->AddLast(track);
 567       }
 568     }
 569   }
 570 }
 571
 572 //_____________________________________________________________________________
 573 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
 574   //-----------------------------------------------------------------
 575   // This function reades track seeds.
 576   //-----------------------------------------------------------------
 577   TDirectory *savedir=gDirectory;
 578
 579   TFile *in=(TFile*)inp;
 580   if (!in->IsOpen()) {
 581      cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
 582      return 1;
 583   }
 584
 585   in->cd();
 586   TTree *seedTree=(TTree*)in->Get("Seeds");
 587   if (!seedTree) {
 588      cerr<<"AliTPCtracker::ReadSeeds(): ";
 589      cerr<<"can't get a tree with track seeds !\n";
 590      return 2;
 591   }
 592   AliTPCtrack *seed=new AliTPCtrack;
 593   seedTree->SetBranchAddress("tracks",&seed);
 594
 595   if (fSeeds==0) fSeeds=new TObjArray(15000);
 596
 597   Int_t n=(Int_t)seedTree->GetEntries();
 598   for (Int_t i=0; i<n; i++) {
 599      seedTree->GetEvent(i);
 600      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
 601   }
 602
 603   delete seed;
 604
 605   delete seedTree; //Thanks to Mariana Bondila
 606
 607   savedir->cd();
 608
 609   return 0;
 610 }
 611
 612 //_____________________________________________________________________________
 613 Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
 614   //-----------------------------------------------------------------
 615   // This is a track finder.
 616   //-----------------------------------------------------------------
 617   TDirectory *savedir=gDirectory;
 618
 619   if (inp) {
 620      TFile *in=(TFile*)inp;
 621      if (!in->IsOpen()) {
 622         cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
 623         return 1;
 624      }
 625   }
 626
 627   if (!out->IsOpen()) {
 628      cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
 629      return 2;
 630   }
 631
 632   out->cd();
 633
 634   char   tname[100];
 635   sprintf(tname,"TreeT_TPC_%d",fEventN);
 636   TTree tracktree(tname,"Tree with TPC tracks");
 637   AliTPCtrack *iotrack=0;
 638   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 639
 640   LoadOuterSectors();
 641
 642   //find track seeds
 643   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 644   Int_t nrows=nlow+nup;
 645   if (fSeeds==0) {
 646      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
 647      MakeSeeds(nup-1, nup-1-gap);
 648      MakeSeeds(nup-1-shift, nup-1-shift-gap);
 649   }
 650   fSeeds->Sort();
 651
 652   //tracking in outer sectors
 653   Int_t nseed=fSeeds->GetEntriesFast();
 654   Int_t i;
 655   for (i=0; i<nseed; i++) {
 656     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 657     if (FollowProlongation(t)) {
 658        UseClusters(&t);
 659        continue;
 660     }
 661     delete fSeeds->RemoveAt(i);
 662   }
 663   //UnloadOuterSectors();
 664
 665   //tracking in inner sectors
 666   LoadInnerSectors();
 667   Int_t found=0;
 668   for (i=0; i<nseed; i++) {
 669     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 670     if (!pt) continue;
 671     Int_t nc=t.GetNumberOfClusters();
 672
 673     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
 674     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 675     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 676     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
 677
 678     alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
 679
 680     if (t.Rotate(alpha)) {
 681        if (FollowProlongation(t)) {
 682           if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 683              t.CookdEdx();
 684              CookLabel(pt,0.1); //For comparison only
 685              iotrack=pt;
 686              tracktree.Fill();
 687              UseClusters(&t,nc);
 688              cerr<<found++<<'\r';
 689           }
 690        }
 691     }
 692     delete fSeeds->RemoveAt(i);
 693   }
 694   UnloadInnerSectors();
 695   UnloadOuterSectors();
 696
 697   tracktree.Write();
 698
 699   cerr<<"Number of found tracks : "<<found<<endl;
 700
 701   savedir->cd();
 702
 703   return 0;
 704 }
 705
 706 //_____________________________________________________________________________
 707 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
 708   //-----------------------------------------------------------------
 709   // This function propagates tracks back through the TPC.
 710   //-----------------------------------------------------------------
 711   fSeeds=new TObjArray(15000);
 712
 713   TFile *in=(TFile*)inp;
 714   TDirectory *savedir=gDirectory;
 715
 716   if (!in->IsOpen()) {
 717      cerr<<"AliTPCtracker::PropagateBack(): ";
 718      cerr<<"file with back propagated ITS tracks is not open !\n";
 719      return 1;
 720   }
 721
 722   if (!out->IsOpen()) {
 723      cerr<<"AliTPCtracker::PropagateBack(): ";
 724      cerr<<"file for back propagated TPC tracks is not open !\n";
 725      return 2;
 726   }
 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   LoadInnerSectors();
 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   UnloadInnerSectors();
 793
 794 //*** Back propagation through outer sectors
 795   LoadOuterSectors();
 796   Int_t found=0;
 797   for (i=0; i<nseed; i++) {
 798      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
 799      if (!ps) continue;
 800      Int_t nc=s.GetNumberOfClusters();
 801      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
 802
 803      Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
 804      if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 805      if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 806      Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
 807
 808      alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
 809      alpha-=s.GetAlpha();
 810
 811      if (s.Rotate(alpha)) {
 812         if (FollowBackProlongation(s,t)) {
 813            if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 814               s.CookdEdx();
 815               s.SetLabel(t.GetLabel());
 816               UseClusters(&s,nc);
 817               otrack=ps;
 818               backTree.Fill();
 819               cerr<<found++<<'\r';
 820               continue;
 821            }
 822         }
 823      }
 824      delete fSeeds->RemoveAt(i);
 825   }
 826   UnloadOuterSectors();
 827
 828   backTree.Write();
 829   savedir->cd();
 830   cerr<<"Number of seeds: "<<nseed<<endl;
 831   cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
 832   cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
 833
 834   delete bckTrack;
 835   delete tpcTrack;
 836
 837   delete bckTree; //Thanks to Mariana Bondila
 838   delete tpcTree; //Thanks to Mariana Bondila
 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   AliTPCClustersRow *clrow=((AliTPCtracker *) this)->fClustersArray.GetRow(sec,row);
 853   return (AliCluster*)(*clrow)[ncl];
 854 }
 855
 856 //__________________________________________________________________________
 857 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
 858   //--------------------------------------------------------------------
 859   //This function "cooks" a track label. If label<0, this track is fake.
 860   //--------------------------------------------------------------------
 861   Int_t noc=t->GetNumberOfClusters();
 862   Int_t *lb=new Int_t[noc];
 863   Int_t *mx=new Int_t[noc];
 864   AliCluster **clusters=new AliCluster*[noc];
 865
 866   Int_t i;
 867   for (i=0; i<noc; i++) {
 868      lb[i]=mx[i]=0;
 869      Int_t index=t->GetClusterIndex(i);
 870      clusters[i]=GetCluster(index);
 871   }
 872
 873   Int_t lab=123456789;
 874   for (i=0; i<noc; i++) {
 875     AliCluster *c=clusters[i];
 876     lab=TMath::Abs(c->GetLabel(0));
 877     Int_t j;
 878     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
 879     lb[j]=lab;
 880     (mx[j])++;
 881   }
 882
 883   Int_t max=0;
 884   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
 885
 886   for (i=0; i<noc; i++) {
 887     AliCluster *c=clusters[i];
 888     if (TMath::Abs(c->GetLabel(1)) == lab ||
 889         TMath::Abs(c->GetLabel(2)) == lab ) max++;
 890   }
 891
 892   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
 893
 894   else {
 895      Int_t tail=Int_t(0.10*noc);
 896      max=0;
 897      for (i=1; i<=tail; i++) {
 898        AliCluster *c=clusters[noc-i];
 899        if (lab == TMath::Abs(c->GetLabel(0)) ||
 900            lab == TMath::Abs(c->GetLabel(1)) ||
 901            lab == TMath::Abs(c->GetLabel(2))) max++;
 902      }
 903      if (max < Int_t(0.5*tail)) lab=-lab;
 904   }
 905
 906   t->SetLabel(lab);
 907
 908   delete[] lb;
 909   delete[] mx;
 910   delete[] clusters;
 911 }
 912
 913 //_________________________________________________________________________
 914 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
 915   //-----------------------------------------------------------------------
 916   // Setup inner sector
 917   //-----------------------------------------------------------------------
 918   if (f==0) {
 919      fAlpha=par->GetInnerAngle();
 920      fAlphaShift=par->GetInnerAngleShift();
 921      fPadPitchWidth=par->GetInnerPadPitchWidth();
 922      fPadPitchLength=par->GetInnerPadPitchLength();
 923
 924      fN=par->GetNRowLow();
 925      fRow=new AliTPCRow[fN];
 926      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
 927   } else {
 928      fAlpha=par->GetOuterAngle();
 929      fAlphaShift=par->GetOuterAngleShift();
 930      fPadPitchWidth=par->GetOuterPadPitchWidth();
 931      f1PadPitchLength=par->GetOuter1PadPitchLength();
 932      f2PadPitchLength=par->GetOuter2PadPitchLength();
 933
 934      fN=par->GetNRowUp();
 935      fRow=new AliTPCRow[fN];
 936      for (Int_t i=0; i<fN; i++){
 937      fRow[i].SetX(par->GetPadRowRadiiUp(i));
 938      }
 939   }
 940 }
 941
 942 //_________________________________________________________________________
 943 void
 944 AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) {
 945   //-----------------------------------------------------------------------
 946   // Insert a cluster into this pad row in accordence with its y-coordinate
 947   //-----------------------------------------------------------------------
 948   if (fN==kMaxClusterPerRow) {
 949     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
 950   }
 951   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
 952   Int_t i=Find(c->GetY());
 953   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
 954   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
 955   fIndex[i]=index; fClusters[i]=c; fN++;
 956 }
 957
 958 //___________________________________________________________________
 959 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
 960   //-----------------------------------------------------------------------
 961   // Return the index of the nearest cluster
 962   //-----------------------------------------------------------------------
 963   if (y <= fClusters[0]->GetY()) return 0;
 964   if (y > fClusters[fN-1]->GetY()) return fN;
 965   Int_t b=0, e=fN-1, m=(b+e)/2;
 966   for (; b<e; m=(b+e)/2) {
 967     if (y > fClusters[m]->GetY()) b=m+1;
 968     else e=m;
 969   }
 970   return m;
 971 }
 972
 973 //_____________________________________________________________________________
 974 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
 975   //-----------------------------------------------------------------
 976   // This funtion calculates dE/dX within the "low" and "up" cuts.
 977   //-----------------------------------------------------------------
 978   Int_t i;
 979   Int_t nc=GetNumberOfClusters();
 980
 981   Int_t swap;//stupid sorting
 982   do {
 983     swap=0;
 984     for (i=0; i<nc-1; i++) {
 985       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
 986       Float_t tmp=fdEdxSample[i];
 987       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
 988       swap++;
 989     }
 990   } while (swap);
 991
 992   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
 993   Float_t dedx=0;
 994   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
 995   dedx /= (nu-nl+1);
 996   SetdEdx(dedx);
 997
 998   //Very rough PID
 999   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1000
1001   if (p<0.6) {
1002     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1003     if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
1004     SetMass(0.93827); return;
1005   }
1006
1007   if (p<1.2) {
1008     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1009     SetMass(0.93827); return;
1010   }
1011
1012   SetMass(0.13957); return;
1013
1014 }
1015
1016
1017