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