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