TPC/AliTPCtracker.cxx

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