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.32  2003/04/10 10:36:54  hristov
  19 Code for unified TPC/TRD tracking (S.Radomski)
  20
  21 Revision 1.31  2003/03/19 17:14:11  hristov
  22 Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
  23
  24 Revision 1.30  2003/02/28 16:13:32  hristov
  25 Typos corrected
  26
  27 Revision 1.29  2003/02/28 15:18:16  hristov
  28 Corrections suggested by J.Chudoba
  29
  30 Revision 1.28  2003/02/27 16:15:52  hristov
  31 Code for inward refitting (S.Radomski)
  32
  33 Revision 1.27  2003/02/25 16:47:58  hristov
  34 allow back propagation for more than 1 event (J.Chudoba)
  35
  36 Revision 1.26  2003/02/19 08:49:46  hristov
  37 Track time measurement (S.Radomski)
  38
  39 Revision 1.25  2003/01/28 16:43:35  hristov
  40 Additional protection: to be discussed with the Root team (M.Ivanov)
  41
  42 Revision 1.24  2002/11/19 16:13:24  hristov
  43 stdlib.h included to declare exit() on HP
  44
  45 Revision 1.23  2002/11/19 11:50:08  hristov
  46 Removing CONTAINERS (Yu.Belikov)
  47
  48 Revision 1.19  2002/07/19 07:31:40  kowal2
  49 Improvement in tracking by J. Belikov
  50
  51 Revision 1.18  2002/05/13 07:33:52  kowal2
  52 Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
  53 in the case of defined region of interests.
  54
  55 Revision 1.17  2002/03/18 17:59:13  kowal2
  56 Chnges in the pad geometry - 3 pad lengths introduced.
  57
  58 Revision 1.16  2001/11/08 16:39:03  hristov
  59 Additional protection (M.Masera)
  60
  61 Revision 1.15  2001/11/08 16:36:33  hristov
  62 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)
  63
  64 Revision 1.14  2001/10/21 19:04:55  hristov
  65 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
  66
  67 Revision 1.13  2001/07/23 12:01:30  hristov
  68 Initialisation added
  69
  70 Revision 1.12  2001/07/20 14:32:44  kowal2
  71 Processing of many events possible now
  72
  73 Revision 1.11  2001/05/23 08:50:10  hristov
  74 Weird inline removed
  75
  76 Revision 1.10  2001/05/16 14:57:25  alibrary
  77 New files for folders and Stack
  78
  79 Revision 1.9  2001/05/11 07:16:56  hristov
  80 Fix needed on Sun and Alpha
  81
  82 Revision 1.8  2001/05/08 15:00:15  hristov
  83 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
  84
  85 Revision 1.5  2000/12/20 07:51:59  kowal2
  86 Changes suggested by Alessandra and Paolo to avoid overlapped
  87 data fields in encapsulated classes.
  88
  89 Revision 1.4  2000/11/02 07:27:16  kowal2
  90 code corrections
  91
  92 Revision 1.2  2000/06/30 12:07:50  kowal2
  93 Updated from the TPC-PreRelease branch
  94
  95 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
  96 Splitted from AliTPCtracking
  97
  98 */
  99
 100 //-------------------------------------------------------
 101 //          Implementation of the TPC tracker
 102 //
 103 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 104 //-------------------------------------------------------
 105 #include <TObjArray.h>
 106 #include <TFile.h>
 107 #include <TTree.h>
 108 #include <Riostream.h>
 109
 110 #include "AliESD.h"
 111
 112 #include "AliTPCtracker.h"
 113 #include "AliTPCcluster.h"
 114 #include "AliTPCParam.h"
 115 #include "AliClusters.h"
 116
 117 #include "TVector2.h"
 118
 119 #include <stdlib.h>
 120
 121 ClassImp(AliTPCtracker)
 122
 123 //_____________________________________________________________________________
 124 AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
 125 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 126 {
 127   //---------------------------------------------------------------------
 128   // The main TPC tracker constructor
 129   //---------------------------------------------------------------------
 130   fInnerSec=new AliTPCSector[fkNIS];
 131   fOuterSec=new AliTPCSector[fkNOS];
 132
 133   Int_t i;
 134   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
 135   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
 136
 137   fParam = (AliTPCParam*) par;
 138   fSeeds=0;
 139
 140   // [SR 17.03.2003]
 141
 142   fBarrelFile = 0;
 143   fBarrelTree = 0;
 144   fBarrelArray = 0;
 145   fBarrelTrack = 0;
 146
 147 }
 148
 149 //_____________________________________________________________________________
 150 AliTPCtracker::~AliTPCtracker() {
 151   //------------------------------------------------------------------
 152   // TPC tracker destructor
 153   //------------------------------------------------------------------
 154   delete[] fInnerSec;
 155   delete[] fOuterSec;
 156   if (fSeeds) {
 157     fSeeds->Delete();
 158     delete fSeeds;
 159   }
 160
 161   // [SR, 01.04.2003]
 162   if (fBarrelFile) {
 163     fBarrelFile->Close();
 164     delete fBarrelFile;
 165   }
 166 }
 167
 168 //_____________________________________________________________________________
 169
 170 void AliTPCtracker::SetBarrelTree(const char *mode) {
 171   //
 172   // Creates a tree for BarrelTracks
 173   // mode = "back" or "refit"
 174   //
 175   // [SR, 01.04.2003]
 176   //
 177
 178   if (!IsStoringBarrel()) return;
 179
 180   TDirectory *sav = gDirectory;
 181   if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
 182
 183   char buff[40];
 184   sprintf(buff, "BarrelTPC_%d_%s", GetEventNumber(), mode);
 185
 186   fBarrelFile->cd();
 187   fBarrelTree = new TTree(buff, "Barrel TPC tracks");
 188
 189   if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", 4);
 190   for(Int_t i=0; i<4; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
 191
 192   fBarrelTree->Branch("tracks", &fBarrelArray);
 193
 194   sav->cd();
 195 }
 196 //_____________________________________________________________________________
 197
 198 void AliTPCtracker::StoreBarrelTrack(AliTPCtrack *ps, Int_t refPlane, Int_t isIn) {
 199   //
 200   // Stores Track at a given reference plane
 201   //
 202   // refPlane: 1-4
 203   // isIn: 1 - backward, 2 - refit
 204   //
 205
 206   if (!IsStoringBarrel()) return;
 207   if (refPlane < 0 || refPlane > 4) return;
 208   if (isIn > 2) return;
 209
 210   static Int_t nClusters;
 211   static Int_t nWrong;
 212   static Double_t chi2;
 213   static Int_t index;
 214
 215   Int_t newClusters, newWrong;
 216   Double_t newChi2;
 217
 218   if ( (refPlane == 1 && isIn == kTrackBack) ||
 219        (refPlane == 4 && isIn == kTrackRefit) ) {
 220
 221     fBarrelArray->Clear();
 222     nClusters = nWrong = 0;
 223     chi2 = 0.0;
 224     index = 0;
 225   }
 226
 227   // propagate
 228   Double_t refX = 0;
 229   if (refPlane == 1) refX = fParam->GetInnerRadiusLow();
 230   if (refPlane == 2) refX = fParam->GetInnerRadiusUp();
 231   if (refPlane == 3) refX = fParam->GetOuterRadiusLow();
 232   if (refPlane == 4) refX = fParam->GetOuterRadiusUp();
 233
 234   ps->PropagateTo(refX);
 235
 236   fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
 237   ps->GetBarrelTrack(fBarrelTrack);
 238
 239   newClusters = ps->GetNumberOfClusters() - nClusters;
 240   newWrong = ps->GetNWrong() - nWrong;
 241   newChi2 = ps->GetChi2() - chi2;
 242
 243   nClusters =  ps->GetNumberOfClusters();
 244   nWrong = ps->GetNWrong();
 245   chi2 = ps->GetChi2();
 246
 247   fBarrelTrack->SetNClusters(newClusters, newChi2);
 248   fBarrelTrack->SetNWrongClusters(newWrong);
 249   fBarrelTrack->SetRefPlane(refPlane, isIn);
 250 }
 251
 252 //_____________________________________________________________________________
 253 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
 254 {
 255   //
 256   // Parametrised error of the cluster reconstruction (pad direction)
 257   //
 258   // Sigma rphi
 259   const Float_t kArphi=0.41818e-2;
 260   const Float_t kBrphi=0.17460e-4;
 261   const Float_t kCrphi=0.30993e-2;
 262   const Float_t kDrphi=0.41061e-3;
 263
 264   pt=TMath::Abs(pt)*1000.;
 265   Double_t x=r/pt;
 266   tgl=TMath::Abs(tgl);
 267   Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
 268   if (s<0.4e-3) s=0.4e-3;
 269   s*=1.3; //Iouri Belikov
 270
 271   return s;
 272 }
 273
 274 //_____________________________________________________________________________
 275 Double_t SigmaZ2(Double_t r, Double_t tgl)
 276 {
 277   //
 278   // Parametrised error of the cluster reconstruction (drift direction)
 279   //
 280   // Sigma z
 281   const Float_t kAz=0.39614e-2;
 282   const Float_t kBz=0.22443e-4;
 283   const Float_t kCz=0.51504e-1;
 284
 285
 286   tgl=TMath::Abs(tgl);
 287   Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
 288   if (s<0.4e-3) s=0.4e-3;
 289   s*=1.3; //Iouri Belikov
 290
 291   return s;
 292 }
 293
 294 //_____________________________________________________________________________
 295 Double_t f1(Double_t x1,Double_t y1,
 296                    Double_t x2,Double_t y2,
 297                    Double_t x3,Double_t y3)
 298 {
 299   //-----------------------------------------------------------------
 300   // Initial approximation of the track curvature
 301   //-----------------------------------------------------------------
 302   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 303   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 304                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 305   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 306                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 307
 308   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 309
 310   return -xr*yr/sqrt(xr*xr+yr*yr);
 311 }
 312
 313
 314 //_____________________________________________________________________________
 315 Double_t f2(Double_t x1,Double_t y1,
 316                    Double_t x2,Double_t y2,
 317                    Double_t x3,Double_t y3)
 318 {
 319   //-----------------------------------------------------------------
 320   // Initial approximation of the track curvature times center of curvature
 321   //-----------------------------------------------------------------
 322   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 323   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 324                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 325   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 326                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 327
 328   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 329
 330   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 331 }
 332
 333 //_____________________________________________________________________________
 334 Double_t f3(Double_t x1,Double_t y1,
 335                    Double_t x2,Double_t y2,
 336                    Double_t z1,Double_t z2)
 337 {
 338   //-----------------------------------------------------------------
 339   // Initial approximation of the tangent of the track dip angle
 340   //-----------------------------------------------------------------
 341   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 342 }
 343
 344 Int_t AliTPCtracker::LoadClusters() {
 345   return LoadClusters(gFile);
 346 }
 347
 348 //_____________________________________________________________________________
 349 Int_t AliTPCtracker::LoadClusters(const TFile *cf) {
 350   //-----------------------------------------------------------------
 351   // This function loads TPC clusters.
 352   //-----------------------------------------------------------------
 353   if (!((TFile*)cf)->IsOpen()) {
 354     cerr<<"AliTPCtracker::LoadClusters : "<<
 355       "file with clusters has not been open !\n";
 356     return 1;
 357   }
 358
 359   Char_t name[100];
 360   sprintf(name,"TreeC_TPC_%d",GetEventNumber());
 361   TTree *cTree=(TTree*)((TFile*)cf)->Get(name);
 362   if (!cTree) {
 363     cerr<<"AliTPCtracker::LoadClusters : "<<
 364       "can't get the tree with TPC clusters !\n";
 365     return 2;
 366   }
 367
 368   TBranch *branch=cTree->GetBranch("Segment");
 369   if (!branch) {
 370     cerr<<"AliTPCtracker::LoadClusters : "<<
 371       "can't get the segment branch !\n";
 372     return 3;
 373   }
 374
 375   AliClusters carray, *addr=&carray;
 376   carray.SetClass("AliTPCcluster");
 377   carray.SetArray(0);
 378   branch->SetAddress(&addr);
 379
 380   Int_t nentr=(Int_t)cTree->GetEntries();
 381
 382   for (Int_t i=0; i<nentr; i++) {
 383       cTree->GetEvent(i);
 384
 385       Int_t ncl=carray.GetArray()->GetEntriesFast();
 386
 387       Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
 388       Int_t id=carray.GetID();
 389       if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
 390          cerr<<"AliTPCtracker::LoadClusters : "<<
 391                "wrong index !\n";
 392          exit(1);
 393       }
 394       Int_t outindex = 2*fkNIS*nir;
 395       if (id<outindex) {
 396          Int_t sec = id/nir;
 397          Int_t row = id - sec*nir;
 398          sec %= fkNIS;
 399          AliTPCRow &padrow=fInnerSec[sec][row];
 400          while (ncl--) {
 401            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
 402            padrow.InsertCluster(c,sec,row);
 403          }
 404       } else {
 405          id -= outindex;
 406          Int_t sec = id/nor;
 407          Int_t row = id - sec*nor;
 408          sec %= fkNOS;
 409          AliTPCRow &padrow=fOuterSec[sec][row];
 410          while (ncl--) {
 411            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
 412            padrow.InsertCluster(c,sec+fkNIS,row);
 413          }
 414       }
 415       carray.GetArray()->Clear();
 416   }
 417   delete cTree;
 418   return 0;
 419 }
 420
 421 //_____________________________________________________________________________
 422 void AliTPCtracker::UnloadClusters() {
 423   //-----------------------------------------------------------------
 424   // This function unloads TPC clusters.
 425   //-----------------------------------------------------------------
 426   Int_t i;
 427   for (i=0; i<fkNIS; i++) {
 428     Int_t nr=fInnerSec->GetNRows();
 429     for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
 430   }
 431   for (i=0; i<fkNOS; i++) {
 432     Int_t nr=fOuterSec->GetNRows();
 433     for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
 434   }
 435 }
 436
 437 //_____________________________________________________________________________
 438 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
 439   //-----------------------------------------------------------------
 440   // This function tries to find a track prolongation.
 441   //-----------------------------------------------------------------
 442   Double_t xt=t.GetX();
 443   const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip :
 444                     Int_t(0.5*fSectors->GetNRows());
 445   Int_t tryAgain=kSKIP;
 446
 447   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
 448   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 449   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 450   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 451
 452   Int_t nrows=fSectors->GetRowNumber(xt)-1;
 453   for (Int_t nr=nrows; nr>=rf; nr--) {
 454     Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
 455     if (!t.PropagateTo(x)) return 0;
 456
 457     AliTPCcluster *cl=0;
 458     UInt_t index=0;
 459     Double_t maxchi2=kMaxCHI2;
 460     const AliTPCRow &krow=fSectors[s][nr];
 461     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
 462     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
 463     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
 464     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 465
 466     if (road>kMaxROAD) {
 467       if (t.GetNumberOfClusters()>4)
 468         cerr<<t.GetNumberOfClusters()
 469         <<"FindProlongation warning: Too broad road !\n";
 470       return 0;
 471     }
 472
 473     //if (fSectors == fInnerSec && (nr == 63 || nr == 0)) {
 474     //  cout << nr << "\t" << krow << endl;
 475     //}
 476
 477     if (krow) {
 478       for (Int_t i=krow.Find(y-road); i<krow; i++) {
 479         AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 480         if (c->GetY() > y+road) break;
 481         if (c->IsUsed()) continue;
 482         if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
 483         Double_t chi2=t.GetPredictedChi2(c);
 484         if (chi2 > maxchi2) continue;
 485         maxchi2=chi2;
 486         cl=c;
 487         index=krow.GetIndex(i);
 488       }
 489     }
 490     if (cl) {
 491       Float_t l=fSectors->GetPadPitchWidth();
 492       Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
 493       t.SetSampledEdx(cl->GetQ()/l*corr,t.GetNumberOfClusters());
 494       if (!t.Update(cl,maxchi2,index)) {
 495          if (!tryAgain--) return 0;
 496       } else tryAgain=kSKIP;
 497     } else {
 498       if (tryAgain==0) break;
 499       if (y > ymax) {
 500          s = (s+1) % fN;
 501          if (!t.Rotate(fSectors->GetAlpha())) return 0;
 502       } else if (y <-ymax) {
 503          s = (s-1+fN) % fN;
 504          if (!t.Rotate(-fSectors->GetAlpha())) return 0;
 505       }
 506       tryAgain--;
 507     }
 508   }
 509
 510   return 1;
 511 }
 512 //_____________________________________________________________________________
 513
 514 Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
 515   //
 516   // This function propagates seed inward TPC using old clusters
 517   // from the track.
 518   //
 519   // Sylwester Radomski, GSI
 520   // 26.02.2003
 521   //
 522
 523   Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
 524   TVector2::Phi_0_2pi(alpha);
 525   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 526
 527   Int_t idx=-1, sec=-1, row=-1;
 528   Int_t nc = seed->GetNumber();
 529
 530   idx=track->GetClusterIndex(nc);
 531   sec=(idx&0xff000000)>>24;
 532   row=(idx&0x00ff0000)>>16;
 533
 534   if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 535   else { if (sec <  fkNIS) row=-1; }
 536
 537   Int_t nr=fSectors->GetNRows();
 538   for (Int_t i=nr-1; i>=0; i--) {
 539
 540     //cout << endl;
 541     //cout << i << "\t" << nc << "\t"<< row << "\t" << seed->Pt() << "\t" << track->Pt() << "\t\t";
 542
 543     Double_t x=fSectors->GetX(i);
 544     if (!seed->PropagateTo(x)) return 0;
 545
 546     // Change sector and rotate seed if necessary
 547
 548     Double_t ymax=fSectors->GetMaxY(i);
 549     Double_t y=seed->GetY();
 550
 551     if (y > ymax) {
 552        s = (s+1) % fN;
 553        if (!seed->Rotate(fSectors->GetAlpha())) return 0;
 554     } else if (y <-ymax) {
 555        s = (s-1+fN) % fN;
 556        if (!seed->Rotate(-fSectors->GetAlpha())) return 0;
 557     }
 558
 559     // try to find a cluster
 560
 561     AliTPCcluster *cl=0;
 562     Int_t index = 0;
 563     Double_t maxchi2 = kMaxCHI2;
 564
 565     if (row==i) {
 566
 567       //cout << row << endl;
 568       // accept already found cluster
 569       AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
 570       Double_t chi2 = seed->GetPredictedChi2(c);
 571
 572       index=idx;
 573       cl=c;
 574       maxchi2=chi2;
 575
 576       //cout << c->GetY() << "\t" << seed->GetY() << "\t" << track->GetY();
 577
 578       if (--nc >= 0 /* track->GetNumberOfClusters()*/ ) {
 579         idx=track->GetClusterIndex(nc);
 580         sec=(idx&0xff000000)>>24;
 581         row=(idx&0x00ff0000)>>16;
 582       }
 583
 584       if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 585       else { if (sec < fkNIS) row=-1; }
 586
 587     }
 588
 589     if (cl) {
 590       Float_t l=fSectors->GetPadPitchWidth();
 591       Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
 592       seed->SetSampledEdx(cl->GetQ()/l*corr,seed->GetNumberOfClusters());
 593       seed->Update(cl,maxchi2,index);
 594     }
 595   }
 596
 597   seed->SetNumber(nc);
 598   return 1;
 599 }
 600
 601 //_____________________________________________________________________________
 602 Int_t AliTPCtracker::FollowBackProlongation
 603 (AliTPCseed& seed, const AliTPCtrack &track) {
 604   //-----------------------------------------------------------------
 605   // This function propagates tracks back through the TPC
 606   //-----------------------------------------------------------------
 607   Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
 608   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 609   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 610   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 611
 612   Int_t idx=-1, sec=-1, row=-1;
 613   //Int_t nc=seed.GetLabel(); //index of the cluster to start with
 614   Int_t nc=seed.GetNumber();
 615
 616   if (nc--) {
 617      idx=track.GetClusterIndex(nc);
 618      sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 619   }
 620   if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 621   else { if (sec <  fkNIS) row=-1; }
 622
 623   Int_t nr=fSectors->GetNRows();
 624   for (Int_t i=0; i<nr; i++) {
 625     Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
 626
 627     //cout << i << "\t" << nc << "\t" << row << "\t" << track.GetNumberOfClusters() << endl;
 628
 629     if (!seed.PropagateTo(x)) return 0;
 630
 631     Double_t y=seed.GetY();
 632     if (y > ymax) {
 633        s = (s+1) % fN;
 634        if (!seed.Rotate(fSectors->GetAlpha())) return 0;
 635     } else if (y <-ymax) {
 636        s = (s-1+fN) % fN;
 637        if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
 638     }
 639
 640     AliTPCcluster *cl=0;
 641     Int_t index=0;
 642     Double_t maxchi2=kMaxCHI2;
 643     Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
 644     Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
 645     Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
 646     Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
 647     if (road>kMaxROAD) {
 648       cerr<<seed.GetNumberOfClusters()
 649           <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n";
 650       return 0;
 651     }
 652
 653     Int_t accepted=seed.GetNumberOfClusters();
 654     if (row==i) {
 655
 656       //if (fSectors == fInnerSec && row == 0)
 657       //cout << "row == " << row << endl;
 658
 659        //try to accept already found cluster
 660        AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
 661        Double_t chi2;
 662        if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) {
 663           index=idx; cl=c; maxchi2=chi2;
 664        } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
 665
 666        if (nc--) {
 667           idx=track.GetClusterIndex(nc);
 668           sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
 669        }
 670        if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
 671        else { if (sec < fkNIS) row=-1; }
 672
 673     }
 674     if (!cl) {
 675        //try to fill the gap
 676        const AliTPCRow &krow=fSectors[s][i];
 677        if (accepted>27)
 678        if (krow) {
 679           for (Int_t i=krow.Find(y-road); i<krow; i++) {
 680             AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
 681             if (c->GetY() > y+road) break;
 682             if (c->IsUsed()) continue;
 683          if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
 684             Double_t chi2=seed.GetPredictedChi2(c);
 685             if (chi2 > maxchi2) continue;
 686             maxchi2=chi2;
 687             cl=c;
 688             index=krow.GetIndex(i);
 689           }
 690        }
 691     }
 692
 693     if (cl) {
 694       Float_t l=fSectors->GetPadPitchWidth();
 695       Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
 696       seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
 697       seed.Update(cl,maxchi2,index);
 698     }
 699
 700   }
 701
 702   //seed.SetLabel(nc);
 703   seed.SetNumber(nc);
 704
 705   return 1;
 706 }
 707
 708 //_____________________________________________________________________________
 709 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
 710   //-----------------------------------------------------------------
 711   // This function creates track seeds.
 712   //-----------------------------------------------------------------
 713   Double_t x[5], c[15];
 714
 715   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
 716   Double_t cs=cos(alpha), sn=sin(alpha);
 717
 718   Double_t x1 =fSectors->GetX(i1);
 719   Double_t xx2=fSectors->GetX(i2);
 720
 721   for (Int_t ns=0; ns<fN; ns++) {
 722     Int_t nl=fSectors[(ns-1+fN)%fN][i2];
 723     Int_t nm=fSectors[ns][i2];
 724     Int_t nu=fSectors[(ns+1)%fN][i2];
 725     const AliTPCRow& kr1=fSectors[ns][i1];
 726     for (Int_t is=0; is < kr1; is++) {
 727       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
 728       for (Int_t js=0; js < nl+nm+nu; js++) {
 729         const AliTPCcluster *kcl;
 730         Double_t x2,   y2,   z2;
 731         Double_t x3=GetX(), y3=GetY(), z3=GetZ();
 732
 733         if (js<nl) {
 734           const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
 735           kcl=kr2[js];
 736           y2=kcl->GetY(); z2=kcl->GetZ();
 737           x2= xx2*cs+y2*sn;
 738           y2=-xx2*sn+y2*cs;
 739         } else
 740           if (js<nl+nm) {
 741             const AliTPCRow& kr2=fSectors[ns][i2];
 742             kcl=kr2[js-nl];
 743             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
 744           } else {
 745             const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
 746             kcl=kr2[js-nl-nm];
 747             y2=kcl->GetY(); z2=kcl->GetZ();
 748             x2=xx2*cs-y2*sn;
 749             y2=xx2*sn+y2*cs;
 750           }
 751
 752         Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
 753         if (TMath::Abs(zz-z2)>5.) continue;
 754
 755         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
 756         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
 757
 758         x[0]=y1;
 759         x[1]=z1;
 760         x[4]=f1(x1,y1,x2,y2,x3,y3);
 761         if (TMath::Abs(x[4]) >= 0.0066) continue;
 762         x[2]=f2(x1,y1,x2,y2,x3,y3);
 763         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
 764         x[3]=f3(x1,y1,x2,y2,z1,z2);
 765         if (TMath::Abs(x[3]) > 1.2) continue;
 766         Double_t a=asin(x[2]);
 767         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
 768         if (TMath::Abs(zv-z3)>10.) continue;
 769
 770         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
 771         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
 772         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
 773         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
 774
 775         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
 776         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
 777         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
 778         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
 779         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
 780         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
 781         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
 782         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
 783         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
 784         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
 785
 786         c[0]=sy1;
 787         c[1]=0.;       c[2]=sz1;
 788         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
 789         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
 790                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
 791         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
 792         c[13]=f30*sy1*f40+f32*sy2*f42;
 793         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
 794
 795         UInt_t index=kr1.GetIndex(is);
 796         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
 797         Float_t l=fSectors->GetPadPitchWidth();
 798         track->SetSampledEdx(kr1[is]->GetQ()/l,0);
 799
 800         Int_t rc=FollowProlongation(*track, i2);
 801         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
 802         else fSeeds->AddLast(track);
 803       }
 804     }
 805   }
 806 }
 807
 808 //_____________________________________________________________________________
 809 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
 810   //-----------------------------------------------------------------
 811   // This function reades track seeds.
 812   //-----------------------------------------------------------------
 813   TDirectory *savedir=gDirectory;
 814
 815   TFile *in=(TFile*)inp;
 816   if (!in->IsOpen()) {
 817      cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
 818      return 1;
 819   }
 820
 821   in->cd();
 822   TTree *seedTree=(TTree*)in->Get("Seeds");
 823   if (!seedTree) {
 824      cerr<<"AliTPCtracker::ReadSeeds(): ";
 825      cerr<<"can't get a tree with track seeds !\n";
 826      return 2;
 827   }
 828   AliTPCtrack *seed=new AliTPCtrack;
 829   seedTree->SetBranchAddress("tracks",&seed);
 830
 831   if (fSeeds==0) fSeeds=new TObjArray(15000);
 832
 833   Int_t n=(Int_t)seedTree->GetEntries();
 834   for (Int_t i=0; i<n; i++) {
 835      seedTree->GetEvent(i);
 836      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
 837   }
 838
 839   delete seed;
 840
 841   delete seedTree; //Thanks to Mariana Bondila
 842
 843   savedir->cd();
 844
 845   return 0;
 846 }
 847
 848 //_____________________________________________________________________________
 849 Int_t AliTPCtracker::Clusters2Tracks(AliESD *event) {
 850   //-----------------------------------------------------------------
 851   // This is a track finder.
 852   // The clusters must be already loaded !
 853   //-----------------------------------------------------------------
 854
 855   //find track seeds
 856   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 857   Int_t nrows=nlow+nup;
 858   if (fSeeds==0) {
 859      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
 860      fSectors=fOuterSec; fN=fkNOS;
 861      fSeeds=new TObjArray(15000);
 862      MakeSeeds(nup-1, nup-1-gap);
 863      MakeSeeds(nup-1-shift, nup-1-shift-gap);
 864   }
 865   fSeeds->Sort();
 866
 867   Int_t nseed=fSeeds->GetEntriesFast();
 868   for (Int_t i=0; i<nseed; i++) {
 869     //tracking in the outer sectors
 870     fSectors=fOuterSec; fN=fkNOS;
 871
 872     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 873     if (!FollowProlongation(t)) {
 874        delete fSeeds->RemoveAt(i);
 875        continue;
 876     }
 877
 878     //tracking in the inner sectors
 879     fSectors=fInnerSec; fN=fkNIS;
 880
 881     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
 882     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 883     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 884     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
 885
 886     alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
 887
 888     if (t.Rotate(alpha)) {
 889       if (FollowProlongation(t)) {
 890         if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 891           t.CookdEdx();
 892           CookLabel(pt,0.1); //For comparison only
 893           pt->PropagateTo(fParam->GetInnerRadiusLow());
 894           AliESDtrack iotrack;
 895           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 896
 897           event->AddTrack(&iotrack);
 898
 899           UseClusters(&t);
 900         }
 901       }
 902     }
 903     delete fSeeds->RemoveAt(i);
 904   }
 905
 906   cerr<<"Number of found tracks : "<<event->GetNumberOfTracks()<<endl;
 907
 908   fSeeds->Clear(); delete fSeeds; fSeeds=0;
 909
 910   return 0;
 911 }
 912
 913 //_____________________________________________________________________________
 914 Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
 915   //-----------------------------------------------------------------
 916   // This is a track finder.
 917   //-----------------------------------------------------------------
 918   TDirectory *savedir=gDirectory;
 919
 920   if (inp) {
 921      TFile *in=(TFile*)inp;
 922      if (!in->IsOpen()) {
 923         cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
 924         return 1;
 925      }
 926   }
 927
 928   if (!out->IsOpen()) {
 929      cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
 930      return 2;
 931   }
 932
 933   LoadClusters();
 934
 935   out->cd();
 936
 937   char   tname[100];
 938   sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
 939   TTree tracktree(tname,"Tree with TPC tracks");
 940   AliTPCtrack *iotrack=0;
 941   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,3);
 942
 943   //find track seeds
 944   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
 945   Int_t nrows=nlow+nup;
 946   if (fSeeds==0) {
 947      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
 948      fSectors=fOuterSec; fN=fkNOS;
 949      fSeeds=new TObjArray(15000);
 950      MakeSeeds(nup-1, nup-1-gap);
 951      MakeSeeds(nup-1-shift, nup-1-shift-gap);
 952   }
 953   fSeeds->Sort();
 954
 955   Int_t found=0;
 956   Int_t nseed=fSeeds->GetEntriesFast();
 957
 958   cout << fInnerSec->GetNRows() << " " << fOuterSec->GetNRows() << endl;
 959
 960   for (Int_t i=0; i<nseed; i++) {
 961     //tracking in the outer sectors
 962     fSectors=fOuterSec; fN=fkNOS;
 963
 964     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
 965     if (!FollowProlongation(t)) {
 966        delete fSeeds->RemoveAt(i);
 967        continue;
 968     }
 969
 970     //tracking in the inner sectors
 971     fSectors=fInnerSec; fN=fkNIS;
 972
 973     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
 974     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 975     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 976     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
 977
 978     alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
 979
 980     if (t.Rotate(alpha)) {
 981       if (FollowProlongation(t)) {
 982         if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
 983           t.CookdEdx();
 984           CookLabel(pt,0.1); //For comparison only
 985           pt->PropagateTo(fParam->GetInnerRadiusLow());
 986           iotrack=pt;
 987           tracktree.Fill();
 988           UseClusters(&t);
 989           found++;
 990 //          cerr<<found<<'\r';
 991         }
 992       }
 993     }
 994     delete fSeeds->RemoveAt(i);
 995   }
 996
 997   tracktree.Write();
 998
 999   cerr<<"Number of found tracks : "<<found<<endl;
1000
1001   savedir->cd();
1002
1003   UnloadClusters();
1004   fSeeds->Clear(); delete fSeeds; fSeeds=0;
1005
1006   return 0;
1007 }
1008 //_____________________________________________________________________________
1009
1010 Int_t AliTPCtracker::RefitInward(TFile *in, TFile *out) {
1011   //
1012   // The function propagates tracks throught TPC inward
1013   // using already associated clusters.
1014   //
1015   // in - file with back propagated tracks
1016   // outs  - file with inward propagation
1017   //
1018   // Sylwester Radomski, GSI
1019   // 26.02.2003
1020   //
1021
1022   if (!in->IsOpen()) {
1023     cout << "Input File not open !\n" << endl;
1024     return 1;
1025   }
1026
1027   if (!out->IsOpen()) {
1028     cout << "Output File not open !\n" << endl;
1029     return 2;
1030   }
1031
1032   TDirectory *savedir = gDirectory;
1033   LoadClusters();
1034
1035   // Connect to input tree
1036   in->cd();
1037   TTree *inputTree = (TTree*)in->Get("seedsTPCtoTRD_0");
1038   TBranch *inBranch = inputTree->GetBranch("tracks");
1039   AliTPCtrack *inTrack = 0;
1040   inBranch->SetAddress(&inTrack);
1041
1042   Int_t nTracks = (Int_t)inputTree->GetEntries();
1043
1044   // Create output tree
1045   out->cd();
1046   AliTPCtrack *outTrack = new AliTPCtrack();
1047   TTree *outputTree = new TTree("tracksTPC_0","Refited TPC tracks");
1048   outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
1049
1050   // [SR, 01.04.2003] - Barrel tracks
1051   if (IsStoringBarrel()) SetBarrelTree("refit");
1052
1053   Int_t nRefited = 0;
1054
1055   for(Int_t t=0; t < nTracks; t++) {
1056
1057     cout << t << "\r";
1058
1059     inputTree->GetEvent(t);
1060     inTrack->ResetCovariance();
1061     AliTPCseed *seed = new AliTPCseed(*inTrack, inTrack->GetAlpha());
1062
1063     /*
1064     //cout << inTrack->GetNumberOfClusters() << endl;
1065     if (inTrack->GetNumberOfClusters() == 158) {
1066       cout << inTrack->GetNumberOfClusters() << endl;
1067       cout << endl;
1068       for(Int_t c=0; c<inTrack->GetNumberOfClusters(); c++) {
1069         Int_t idx = inTrack->GetClusterIndex(c);
1070         Int_t row=(idx&0x00ff0000)>>16;
1071         cout << c << " " << row << endl;
1072       }
1073     }
1074     */
1075
1076     if (IsStoringBarrel()) StoreBarrelTrack(seed, 4, 2);
1077
1078     seed->SetNumber(inTrack->GetNumberOfClusters()-1);
1079     fSectors = fOuterSec;
1080     Int_t res = FollowRefitInward(seed, inTrack);
1081     UseClusters(seed);
1082     Int_t nc = seed->GetNumberOfClusters();
1083
1084     if (IsStoringBarrel()) StoreBarrelTrack(seed, 3, 2);
1085     if (IsStoringBarrel()) StoreBarrelTrack(seed, 2, 2);
1086
1087     fSectors = fInnerSec;
1088     res = FollowRefitInward(seed, inTrack);
1089     UseClusters(seed, nc);
1090
1091     if (IsStoringBarrel()) StoreBarrelTrack(seed, 1, 2);
1092
1093     if (res) {
1094       seed->PropagateTo(fParam->GetInnerRadiusLow());
1095       outTrack = (AliTPCtrack*)seed;
1096       outTrack->SetLabel(inTrack->GetLabel());
1097       outputTree->Fill();
1098       nRefited++;
1099     }
1100
1101     if (IsStoringBarrel()) fBarrelTree->Fill();
1102     delete seed;
1103   }
1104
1105   cout << "Refitted " << nRefited << " tracks." << endl;
1106
1107   out->cd();
1108   outputTree->Write();
1109
1110   // [SR, 01.04.2003]
1111   if (IsStoringBarrel()) {
1112     fBarrelFile->cd();
1113     fBarrelTree->Write();
1114     fBarrelFile->Flush();
1115
1116     delete fBarrelArray;
1117     delete fBarrelTree;
1118   }
1119
1120   if (inputTree) delete inputTree;
1121   if (outputTree) delete outputTree;
1122
1123   savedir->cd();
1124   UnloadClusters();
1125
1126   return 0;
1127 }
1128
1129 //_____________________________________________________________________________
1130 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
1131   //-----------------------------------------------------------------
1132   // This function propagates tracks back through the TPC.
1133   //-----------------------------------------------------------------
1134   return PropagateBack(inp, NULL, out);
1135 }
1136
1137 Int_t AliTPCtracker::PropagateBack(AliESD *event) {
1138   //-----------------------------------------------------------------
1139   // This function propagates tracks back through the TPC.
1140   // The clusters must be already loaded !
1141   //-----------------------------------------------------------------
1142   Int_t nentr=event->GetNumberOfTracks();
1143   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
1144
1145   Int_t ntrk=0;
1146   for (Int_t i=0; i<nentr; i++) {
1147     AliESDtrack *esd=event->GetTrack(i);
1148     ULong_t status=esd->GetStatus();
1149
1150     if ( (status & AliESDtrack::kTPCin ) == 0 ) continue;
1151     if ( (status & AliESDtrack::kTPCout) != 0 ) continue;
1152
1153     const AliTPCtrack t(*esd);
1154     AliTPCseed s(t,t.GetAlpha());
1155
1156     if (status==AliESDtrack::kTPCin) s.ResetCovariance();
1157     else if ( (status & AliESDtrack::kITSout) == 0 ) continue;
1158
1159     Int_t nc=t.GetNumberOfClusters();
1160     s.SetLabel(nc-1); //set number of the cluster to start with
1161
1162     //inner sectors
1163     fSectors=fInnerSec; fN=fkNIS;
1164
1165     Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1166     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1167     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1168     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1169     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1170     alpha-=s.GetAlpha();
1171
1172     if (!s.Rotate(alpha)) continue;
1173     if (!FollowBackProlongation(s,t)) continue;
1174
1175
1176     //outer sectors
1177     fSectors=fOuterSec; fN=fkNOS;
1178
1179     alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1180     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1181     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1182     ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1183
1184     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1185     alpha-=s.GetAlpha();
1186
1187     if (!s.Rotate(alpha)) continue;
1188     if (!FollowBackProlongation(s,t)) continue;
1189     {
1190     Int_t nrows=fOuterSec->GetNRows()+fInnerSec->GetNRows();
1191     if (s.GetNumberOfClusters() < Int_t(0.4*nrows)) continue;
1192     }
1193     s.PropagateTo(fParam->GetOuterRadiusUp());
1194     s.CookdEdx();
1195     CookLabel(&s,0.1); //For comparison only
1196     UseClusters(&s);
1197     esd->UpdateTrackParams(&s,AliESDtrack::kTPCout);
1198     ntrk++;
1199   }
1200   cerr<<"Number of back propagated tracks: "<<ntrk<<endl;
1201
1202   return 0;
1203 }
1204
1205 //_____________________________________________________________________________
1206 Int_t AliTPCtracker::PropagateBack(const TFile *inp, const TFile *inp2, TFile *out) {
1207   //-----------------------------------------------------------------
1208   // This function propagates tracks back through the TPC.
1209   //-----------------------------------------------------------------
1210   fSeeds=new TObjArray(15000);
1211
1212   TFile *in=(TFile*)inp;
1213   TFile *in2=(TFile*)inp2;
1214   TDirectory *savedir=gDirectory;
1215
1216   if (!in->IsOpen()) {
1217      cerr<<"AliTPCtracker::PropagateBack(): ";
1218      cerr<<"file with TPC (or back propagated ITS) tracks is not open !\n";
1219      return 1;
1220   }
1221
1222   if (!out->IsOpen()) {
1223      cerr<<"AliTPCtracker::PropagateBack(): ";
1224      cerr<<"file for back propagated TPC tracks is not open !\n";
1225      return 2;
1226   }
1227
1228   LoadClusters();
1229
1230   in->cd();
1231   char tName[100];
1232   sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
1233   TTree *bckTree=(TTree*)in->Get(tName);
1234   if (!bckTree && inp2) bckTree=(TTree*)in2->Get(tName);
1235   if (!bckTree) {
1236      cerr<<"AliTPCtracker::PropagateBack() ";
1237      cerr<<"can't get a tree with back propagated ITS tracks !\n";
1238      //return 3;
1239   }
1240   AliTPCtrack *bckTrack=new AliTPCtrack;
1241   if (bckTree) bckTree->SetBranchAddress("tracks",&bckTrack);
1242
1243   sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
1244   TTree *tpcTree=(TTree*)in->Get(tName);
1245   if (!tpcTree) {
1246      cerr<<"AliTPCtracker::PropagateBack() ";
1247      cerr<<"can't get a tree with TPC tracks !\n";
1248      return 4;
1249   }
1250   AliTPCtrack *tpcTrack=new AliTPCtrack;
1251   tpcTree->SetBranchAddress("tracks",&tpcTrack);
1252
1253   // [SR, 01.04.2003] - Barrel tracks
1254   if (IsStoringBarrel()) SetBarrelTree("back");
1255
1256   //*** Prepare an array of tracks to be back propagated
1257   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1258   Int_t nrows=nlow+nup;
1259
1260   //
1261   // Match ITS tracks with old TPC tracks
1262   // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
1263   //
1264   // the algorithm is linear and uses LUT for sorting
1265   // cost of the algorithm: 2 * number of tracks
1266   //
1267
1268   TObjArray tracks(15000);
1269   Int_t i=0;
1270   Int_t tpcN= (Int_t)tpcTree->GetEntries();
1271   Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
1272
1273   // look up table
1274   const Int_t nLab = 200000; // limit on number of primaries (arbitrary)
1275   Int_t lut[nLab];
1276   for(Int_t i=0; i<nLab; i++) lut[i] = -1;
1277
1278   if (bckTree) {
1279     for(Int_t i=0; i<bckN; i++) {
1280       bckTree->GetEvent(i);
1281       Int_t lab = TMath::Abs(bckTrack->GetLabel());
1282       if (lab < nLab) lut[lab] = i;
1283       else {
1284         cerr << "AliTPCtracker: The size of the LUT is too small\n";
1285       }
1286     }
1287   }
1288
1289   for (Int_t i=0; i<tpcN; i++) {
1290     tpcTree->GetEvent(i);
1291     Double_t alpha=tpcTrack->GetAlpha();
1292
1293     if (!bckTree) {
1294
1295       // No ITS - use TPC track only
1296
1297       tpcTrack->ResetCovariance();
1298       AliTPCseed * seed = new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha());
1299
1300       fSeeds->AddLast(seed);
1301       tracks.AddLast(new AliTPCtrack(*tpcTrack));
1302
1303     } else {
1304
1305       // with ITS
1306       // discard not prolongated tracks (to be discussed)
1307
1308       Int_t lab = TMath::Abs(tpcTrack->GetLabel());
1309       if (lab > nLab || lut[lab] < 0) continue;
1310
1311       bckTree->GetEvent(lut[lab]);
1312       bckTrack->Rotate(alpha-bckTrack->GetAlpha());
1313
1314       fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1315       tracks.AddLast(new AliTPCtrack(*tpcTrack));
1316     }
1317   }
1318
1319   // end of matching
1320
1321   /*
1322   TObjArray tracks(15000);
1323   Int_t i=0, j=0;
1324   Int_t tpcN=(Int_t)tpcTree->GetEntries();
1325   Int_t bckN=(Int_t)bckTree->GetEntries();
1326   if (j<bckN) bckTree->GetEvent(j++);
1327   for (i=0; i<tpcN; i++) {
1328      tpcTree->GetEvent(i);
1329      Double_t alpha=tpcTrack->GetAlpha();
1330      if (j<bckN &&
1331      TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
1332         if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
1333         fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
1334         bckTree->GetEvent(j++);
1335      } else {
1336         tpcTrack->ResetCovariance();
1337         fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
1338      }
1339      tracks.AddLast(new AliTPCtrack(*tpcTrack));
1340   }
1341   */
1342
1343   out->cd();
1344
1345   // tree name seedsTPCtoTRD as expected by TRD and as
1346   // discussed and decided in Strasbourg (May 2002)
1347   // [SR, GSI, 18.02.2003]
1348
1349   sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
1350   TTree backTree(tName,"Tree with back propagated TPC tracks");
1351
1352   AliTPCtrack *otrack=0;
1353   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
1354
1355   //
1356   Int_t nRefPlane;
1357   Int_t found=0;
1358   Int_t nseed=fSeeds->GetEntriesFast();
1359
1360   // loop changed [SR, 01.04.2003]
1361
1362   for (i=0; i<nseed; i++) {
1363
1364     nRefPlane = 1;
1365     if (IsStoringBarrel()) fBarrelArray->Clear();
1366
1367      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
1368      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
1369
1370     ps->ResetNWrong();
1371     ps->ResetNRotation();
1372
1373     if (IsStoringBarrel()) StoreBarrelTrack(ps, 1, 1);
1374
1375     // Load outer sectors
1376     fSectors=fInnerSec;
1377     fN=fkNIS;
1378
1379      Int_t nc=t.GetNumberOfClusters();
1380     //s.SetLabel(nc-1); //set number of the cluster to start with
1381     s.SetNumber(nc);
1382
1383     if (FollowBackProlongation(s,t)) UseClusters(&s);
1384     else {
1385       fSeeds->RemoveAt(i);
1386         continue;
1387      }
1388
1389     if (IsStoringBarrel()) StoreBarrelTrack(ps, 2, 1);
1390     if (IsStoringBarrel()) StoreBarrelTrack(ps, 3, 1);
1391
1392     // Load outer sectors
1393     fSectors=fOuterSec;
1394     fN=fkNOS;
1395
1396     nc=s.GetNumberOfClusters();
1397
1398      Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
1399      if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1400      if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1401      Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
1402
1403      alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
1404      alpha-=s.GetAlpha();
1405
1406      if (s.Rotate(alpha)) {
1407        if (FollowBackProlongation(s,t)) {
1408          if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
1409            s.CookdEdx();
1410            s.SetLabel(t.GetLabel());
1411            UseClusters(&s,nc);
1412
1413            // Propagate to outer reference plane for comparison reasons
1414            // reason for keeping fParam object [SR, GSI, 18.02.2003]
1415
1416            ps->PropagateTo(fParam->GetOuterRadiusUp());
1417            otrack=ps;
1418            backTree.Fill();
1419
1420           if (IsStoringBarrel()) StoreBarrelTrack(ps, 4, 1);
1421           cerr<<found++<<'\r';
1422          }
1423        }
1424      }
1425
1426     if (IsStoringBarrel()) fBarrelTree->Fill();
1427      delete fSeeds->RemoveAt(i);
1428   }
1429
1430
1431   out->cd();
1432   backTree.Write();
1433
1434   // [SR, 01.04.2003]
1435   if (IsStoringBarrel()) {
1436     fBarrelFile->cd();
1437     fBarrelTree->Write();
1438     fBarrelFile->Flush();
1439
1440     delete fBarrelArray;
1441     delete fBarrelTree;
1442   }
1443
1444   savedir->cd();
1445   cerr<<"Number of seeds: "<<nseed<<endl;
1446   cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
1447   cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
1448
1449   delete bckTrack;
1450   delete tpcTrack;
1451
1452   if (bckTree) delete bckTree; //Thanks to Mariana Bondila
1453   delete tpcTree; //Thanks to Mariana Bondila
1454
1455   UnloadClusters();
1456
1457   return 0;
1458 }
1459
1460 //_________________________________________________________________________
1461 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
1462   //--------------------------------------------------------------------
1463   //       Return pointer to a given cluster
1464   //--------------------------------------------------------------------
1465   Int_t sec=(index&0xff000000)>>24;
1466   Int_t row=(index&0x00ff0000)>>16;
1467   Int_t ncl=(index&0x0000ffff)>>00;
1468
1469   const AliTPCcluster *cl=0;
1470
1471   if (sec<fkNIS) {
1472     cl=fInnerSec[sec][row].GetUnsortedCluster(ncl);
1473   } else {
1474     sec-=fkNIS;
1475     cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
1476   }
1477
1478   return (AliCluster*)cl;
1479 }
1480
1481 //__________________________________________________________________________
1482 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1483   //--------------------------------------------------------------------
1484   //This function "cooks" a track label. If label<0, this track is fake.
1485   //--------------------------------------------------------------------
1486   Int_t noc=t->GetNumberOfClusters();
1487   Int_t *lb=new Int_t[noc];
1488   Int_t *mx=new Int_t[noc];
1489   AliCluster **clusters=new AliCluster*[noc];
1490
1491   Int_t i;
1492   for (i=0; i<noc; i++) {
1493      lb[i]=mx[i]=0;
1494      Int_t index=t->GetClusterIndex(i);
1495      clusters[i]=GetCluster(index);
1496   }
1497
1498   Int_t lab=123456789;
1499   for (i=0; i<noc; i++) {
1500     AliCluster *c=clusters[i];
1501     lab=TMath::Abs(c->GetLabel(0));
1502     Int_t j;
1503     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1504     lb[j]=lab;
1505     (mx[j])++;
1506   }
1507
1508   Int_t max=0;
1509   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1510
1511   for (i=0; i<noc; i++) {
1512     AliCluster *c=clusters[i];
1513     if (TMath::Abs(c->GetLabel(1)) == lab ||
1514         TMath::Abs(c->GetLabel(2)) == lab ) max++;
1515   }
1516
1517   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1518
1519   else {
1520      Int_t tail=Int_t(0.10*noc);
1521      max=0;
1522      for (i=1; i<=tail; i++) {
1523        AliCluster *c=clusters[noc-i];
1524        if (lab == TMath::Abs(c->GetLabel(0)) ||
1525            lab == TMath::Abs(c->GetLabel(1)) ||
1526            lab == TMath::Abs(c->GetLabel(2))) max++;
1527      }
1528      if (max < Int_t(0.5*tail)) lab=-lab;
1529   }
1530
1531   t->SetLabel(lab);
1532
1533   delete[] lb;
1534   delete[] mx;
1535   delete[] clusters;
1536 }
1537
1538 //_________________________________________________________________________
1539 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1540   //-----------------------------------------------------------------------
1541   // Setup inner sector
1542   //-----------------------------------------------------------------------
1543   if (f==0) {
1544      fAlpha=par->GetInnerAngle();
1545      fAlphaShift=par->GetInnerAngleShift();
1546      fPadPitchWidth=par->GetInnerPadPitchWidth();
1547      f1PadPitchLength=par->GetInnerPadPitchLength();
1548      f2PadPitchLength=f1PadPitchLength;
1549
1550      fN=par->GetNRowLow();
1551      fRow=new AliTPCRow[fN];
1552      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
1553   } else {
1554      fAlpha=par->GetOuterAngle();
1555      fAlphaShift=par->GetOuterAngleShift();
1556      fPadPitchWidth=par->GetOuterPadPitchWidth();
1557      f1PadPitchLength=par->GetOuter1PadPitchLength();
1558      f2PadPitchLength=par->GetOuter2PadPitchLength();
1559
1560      fN=par->GetNRowUp();
1561      fRow=new AliTPCRow[fN];
1562      for (Int_t i=0; i<fN; i++){
1563      fRow[i].SetX(par->GetPadRowRadiiUp(i));
1564      }
1565   }
1566 }
1567
1568 //_________________________________________________________________________
1569 void AliTPCtracker::
1570 AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
1571   //-----------------------------------------------------------------------
1572   // Insert a cluster into this pad row in accordence with its y-coordinate
1573   //-----------------------------------------------------------------------
1574   if (fN==kMaxClusterPerRow) {
1575     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1576   }
1577
1578   Int_t index=(((sec<<8)+row)<<16)+fN;
1579
1580   if (fN==0) {
1581      fIndex[0]=index;
1582      fClusterArray[0]=*c; fClusters[fN++]=fClusterArray;
1583      return;
1584   }
1585
1586   if (fN==fSize) {
1587      Int_t size=fSize*2;
1588      AliTPCcluster *buff=new AliTPCcluster[size];
1589      memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
1590      for (Int_t i=0; i<fN; i++)
1591         fClusters[i]=buff+(fClusters[i]-fClusterArray);
1592      delete[] fClusterArray;
1593      fClusterArray=buff;
1594      fSize=size;
1595   }
1596
1597   Int_t i=Find(c->GetY());
1598   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
1599   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
1600   fIndex[i]=index;
1601   fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
1602 }
1603
1604 //___________________________________________________________________
1605 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
1606   //-----------------------------------------------------------------------
1607   // Return the index of the nearest cluster
1608   //-----------------------------------------------------------------------
1609   if(fN<=0) return 0;
1610   if (y <= fClusters[0]->GetY()) return 0;
1611   if (y > fClusters[fN-1]->GetY()) return fN;
1612   Int_t b=0, e=fN-1, m=(b+e)/2;
1613   for (; b<e; m=(b+e)/2) {
1614     if (y > fClusters[m]->GetY()) b=m+1;
1615     else e=m;
1616   }
1617   return m;
1618 }
1619
1620 //_____________________________________________________________________________
1621 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1622   //-----------------------------------------------------------------
1623   // This funtion calculates dE/dX within the "low" and "up" cuts.
1624   //-----------------------------------------------------------------
1625   Int_t i;
1626   Int_t nc=GetNumberOfClusters();
1627
1628   Int_t swap;//stupid sorting
1629   do {
1630     swap=0;
1631     for (i=0; i<nc-1; i++) {
1632       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1633       Float_t tmp=fdEdxSample[i];
1634       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1635       swap++;
1636     }
1637   } while (swap);
1638
1639   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
1640   Float_t dedx=0;
1641   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
1642   dedx /= (nu-nl+1);
1643   SetdEdx(dedx);
1644
1645   //Very rough PID
1646   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1647
1648   Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
1649   if (p<0.6) {
1650     if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1651     if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
1652     SetMass(0.93827); return;
1653   }
1654
1655   if (p<1.2) {
1656     if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
1657     SetMass(0.93827); return;
1658   }
1659
1660   SetMass(0.13957); return;
1661
1662 }
1663
1664
1665