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