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