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