TRD/AliTRDtracker.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.9  2001/02/14 18:22:26  cblume
  19 Change in the geometry of the padplane
  20
  21 Revision 1.8  2000/12/20 13:00:44  cblume
  22 Modifications for the HP-compiler
  23
  24 Revision 1.7  2000/12/08 16:07:02  cblume
  25 Update of the tracking by Sergei
  26
  27 Revision 1.6  2000/11/30 17:38:08  cblume
  28 Changes to get in line with new STEER and EVGEN
  29
  30 Revision 1.5  2000/11/14 14:40:27  cblume
  31 Correction for the Sun compiler (kTRUE and kFALSE)
  32
  33 Revision 1.4  2000/11/10 14:57:52  cblume
  34 Changes in the geometry constants for the DEC compiler
  35
  36 Revision 1.3  2000/10/15 23:40:01  cblume
  37 Remove AliTRDconst
  38
  39 Revision 1.2  2000/10/06 16:49:46  cblume
  40 Made Getters const
  41
  42 Revision 1.1.2.2  2000/10/04 16:34:58  cblume
  43 Replace include files by forward declarations
  44
  45 Revision 1.1.2.1  2000/09/22 14:47:52  cblume
  46 Add the tracking code
  47
  48 */
  49
  50 #include <iostream.h>
  51
  52 #include <TFile.h>
  53 #include <TROOT.h>
  54 #include <TBranch.h>
  55 #include <TTree.h>
  56
  57 #include "AliRun.h"
  58 #include "AliTRD.h"
  59 #include "AliTRDgeometry.h"
  60 #include "AliTRDrecPoint.h"
  61 #include "AliTRDcluster.h"
  62 #include "AliTRDtrack.h"
  63 #include "AliTRDtrackingSector.h"
  64 #include "AliTRDtimeBin.h"
  65
  66 #include "AliTRDtracker.h"
  67
  68 ClassImp(AliTRDtracker)
  69
  70   const  Int_t     AliTRDtracker::fSeedGap            = 35;
  71   const  Int_t     AliTRDtracker::fSeedStep           = 5;
  72
  73
  74   const  Float_t   AliTRDtracker::fMinClustersInTrack = 0.5;
  75   const  Float_t   AliTRDtracker::fMinClustersInSeed  = 0.5;
  76   const  Float_t   AliTRDtracker::fSeedDepth          = 0.5;
  77   const  Float_t   AliTRDtracker::fSkipDepth          = 0.2;
  78   const  Float_t   AliTRDtracker::fMaxSeedDeltaZ      = 30.;
  79   const  Float_t   AliTRDtracker::fMaxSeedC           = 0.01;
  80   const  Float_t   AliTRDtracker::fMaxSeedTan         = 1.2;
  81   const  Float_t   AliTRDtracker::fMaxSeedVertexZ     = 200.;
  82   const  Float_t   AliTRDtracker::fLabelFraction      = 0.5;
  83   const  Float_t   AliTRDtracker::fWideRoad           = 30.;
  84
  85   const  Double_t  AliTRDtracker::fMaxChi2            = 12.;
  86   const  Double_t  AliTRDtracker::fSeedErrorSY        = 0.1;
  87   const  Double_t  AliTRDtracker::fSeedErrorSY3       = 2.5;
  88   const  Double_t  AliTRDtracker::fSeedErrorSZ        = 0.1;
  89
  90 //____________________________________________________________________
  91 AliTRDtracker::AliTRDtracker()
  92 {
  93   //
  94   // Default constructor
  95   //
  96
  97   fEvent     = 0;
  98   fGeom      = NULL;
  99
 100   fNclusters = 0;
 101   fClusters  = NULL;
 102   fNseeds    = 0;
 103   fSeeds     = NULL;
 104   fNtracks   = 0;
 105   fTracks    = NULL;
 106
 107 }
 108
 109 //____________________________________________________________________
 110 AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
 111                   :TNamed(name, title)
 112 {
 113   fEvent     = 0;
 114   fGeom      = NULL;
 115
 116   fNclusters = 0;
 117   fClusters  = new TObjArray(2000);
 118   fNseeds    = 0;
 119   fSeeds     = new TObjArray(20000);
 120   fNtracks   = 0;
 121   fTracks    = new TObjArray(10000);
 122
 123 }
 124
 125 //___________________________________________________________________
 126 AliTRDtracker::~AliTRDtracker()
 127 {
 128   delete fClusters;
 129   delete fTracks;
 130   delete fSeeds;
 131   delete fGeom;
 132 }
 133
 134 //___________________________________________________________________
 135 void AliTRDtracker::Clusters2Tracks()
 136 {
 137   Int_t inner, outer;
 138   Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
 139   Int_t nSteps = (Int_t) (fTotalNofTimeBins * fSeedDepth)/fSeedStep;
 140
 141   for(Int_t i=0; i<nSteps; i++) {
 142     printf("step %d out of %d \n", i+1, nSteps);
 143     outer=fTotalNofTimeBins-1-i*fSeedStep; inner=outer-fSeedGap;
 144     MakeSeeds(inner,outer);
 145     FindTracks();
 146   }
 147 }
 148
 149 //_____________________________________________________________________
 150 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
 151 {
 152   // Parametrised "expected" error of the cluster reconstruction in Y
 153
 154   Double_t s = 0.2;
 155   return s;
 156 }
 157
 158 //_____________________________________________________________________
 159 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
 160 {
 161   // Parametrised "expected" error of the cluster reconstruction in Z
 162
 163   // Take an example pad size for the time being, check back
 164   // again with Sergei
 165   Double_t s, pad = fGeom->GetRowPadSize(0,2,0);
 166   s = pad * pad /12.;
 167   return s;
 168 }
 169
 170 //_____________________________________________________________________
 171 inline Double_t f1trd(Double_t x1,Double_t y1,
 172                       Double_t x2,Double_t y2,
 173                       Double_t x3,Double_t y3)
 174 {
 175   // Initial approximation of the track curvature
 176   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 177
 178   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 179   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 180                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 181   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 182                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 183
 184   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 185
 186   return -xr*yr/sqrt(xr*xr+yr*yr);
 187 }
 188
 189 //_____________________________________________________________________
 190 inline Double_t f2trd(Double_t x1,Double_t y1,
 191                       Double_t x2,Double_t y2,
 192                       Double_t x3,Double_t y3)
 193 {
 194   // Initial approximation of the track curvature times center of curvature
 195   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 196
 197   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 198   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 199                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 200   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 201                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 202
 203   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 204
 205   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 206 }
 207
 208 //_____________________________________________________________________
 209 inline Double_t f3trd(Double_t x1,Double_t y1,
 210                       Double_t x2,Double_t y2,
 211                       Double_t z1,Double_t z2)
 212 {
 213   // Initial approximation of the tangent of the track dip angle
 214   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 215
 216   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 217 }
 218
 219
 220 //___________________________________________________________________
 221
 222 Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
 223                             Int_t s, Int_t rf)
 224 {
 225   // Starting from current position on track=t this function tries
 226   // to extrapolate the track up to timeBin=rf and to confirm prolongation
 227   // if a close cluster is found. *sec is a pointer to allocated
 228   // array of sectors, in which the initial sector has index=s.
 229
 230   const Int_t TIME_BINS_TO_SKIP=Int_t(fSkipDepth*sec->GetNtimeBins());
 231   Int_t try_again=TIME_BINS_TO_SKIP;
 232
 233   Double_t alpha=AliTRDgeometry::GetAlpha();
 234
 235   Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
 236
 237   for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
 238     Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
 239     if (!t.PropagateTo(x)) {
 240       cerr<<"Can't propagate to x = "<<x<<endl;
 241       return 0;
 242     }
 243
 244     AliTRDcluster *cl=0;
 245     UInt_t index=0;
 246
 247     Double_t max_chi2=fMaxChi2;
 248
 249     AliTRDtimeBin& time_bin=sec[s][nr];
 250     Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
 251     Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
 252     Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 253
 254     if (road>fWideRoad) {
 255       if (t.GetNclusters() > 4) {
 256         cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
 257       }
 258       return 0;
 259     }
 260
 261     if (time_bin) {
 262       for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
 263         AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
 264         if (c->GetY() > y+road) break;
 265         if (c->IsUsed() > 0) continue;
 266
 267         if((c->GetZ()-z)*(c->GetZ()-z) > 25. + sz2) continue;
 268
 269         Double_t chi2=t.GetPredictedChi2(c);
 270
 271         if (chi2 > max_chi2) continue;
 272         max_chi2=chi2;
 273         cl=c;
 274         index=time_bin.GetIndex(i);
 275       }
 276     }
 277     if (cl) {
 278
 279       //      Float_t l=sec->GetPitch();
 280       //      t.SetSampledEdx(cl->fQ/l,Int_t(t));
 281
 282       t.Update(cl,max_chi2,index);
 283
 284       try_again=TIME_BINS_TO_SKIP;
 285     } else {
 286       if (try_again==0) break;
 287       if (y > ymax) {
 288         cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
 289          s = (s+1) % ns;
 290          if (!t.Rotate(alpha)) {
 291            cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
 292            return 0;
 293          }
 294       } else if (y <-ymax) {
 295         cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
 296          s = (s-1+ns) % ns;
 297          if (!t.Rotate(-alpha)) {
 298            cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
 299            return 0;
 300          }
 301       }
 302       try_again--;
 303     }
 304   }
 305
 306   return 1;
 307 }
 308
 309
 310
 311 //_____________________________________________________________________________
 312 void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
 313 {
 314   // Opens a ROOT-file with TRD-clusters and reads in the cluster-tree
 315
 316   ReadClusters(fClusters, clusterfile);
 317
 318   // get geometry from the file with hits
 319
 320   TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile);
 321   if (!fInputFile) {
 322     printf("AliTRDtracker::Open -- ");
 323     printf("Open the ALIROOT-file %s.\n",hitfile);
 324     fInputFile = new TFile(hitfile);
 325   }
 326   else {
 327     printf("AliTRDtracker::Open -- ");
 328     printf("%s is already open.\n",hitfile);
 329   }
 330
 331   // Get AliRun object from file or create it if not on file
 332
 333   gAlice = (AliRun*) fInputFile->Get("gAlice");
 334   if (gAlice) {
 335     printf("AliTRDtracker::GetEvent -- ");
 336     printf("AliRun object found on file.\n");
 337   }
 338   else {
 339     printf("AliTRDtracker::GetEvent -- ");
 340     printf("Could not find AliRun object.\n");
 341   }
 342
 343   // Import the Trees for the event nEvent in the file
 344   Int_t nparticles = gAlice->GetEvent(fEvent);
 345   cerr<<"nparticles = "<<nparticles<<endl;
 346
 347   if (nparticles <= 0) {
 348     printf("AliTRDtracker::GetEvent -- ");
 349     printf("No entries in the trees for event %d.\n",fEvent);
 350   }
 351
 352   AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
 353   fGeom = TRD->GetGeometry();
 354
 355 }
 356
 357
 358 //_____________________________________________________________________________
 359 void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
 360 {
 361   // Fills clusters into TRD tracking_sectors
 362   // Note that the numbering scheme for the TRD tracking_sectors
 363   // differs from that of TRD sectors
 364
 365   for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp();
 366
 367   //  Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
 368
 369   cerr<<"MakeSeeds: sorting clusters"<<endl;
 370
 371   Int_t ncl=fClusters->GetEntriesFast();
 372   UInt_t index;
 373   while (ncl--) {
 374     printf("\r %d left",ncl);
 375     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
 376     Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
 377     Int_t sector=fGeom->GetSector(detector);
 378
 379     Int_t tracking_sector=sector;
 380     if(sector > 0) tracking_sector = AliTRDgeometry::kNsect - sector;
 381
 382     Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin);
 383     index=ncl;
 384     sec[tracking_sector][tb].InsertCluster(c,index);
 385   }
 386   printf("\r\n");
 387 }
 388
 389
 390 //_____________________________________________________________________________
 391 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
 392 {
 393   // Creates track seeds using clusters in timeBins=i1,i2
 394
 395   Int_t i2 = inner, i1 = outer;
 396
 397   if (!fClusters) return;
 398
 399   AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
 400   SetUpSectors(fTrSec);
 401
 402   // find seeds
 403
 404   Double_t x[5], c[15];
 405   Int_t max_sec=AliTRDgeometry::kNsect;
 406
 407   Double_t alpha=AliTRDgeometry::GetAlpha();
 408   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
 409   Double_t cs=cos(alpha), sn=sin(alpha);
 410
 411   Double_t x1 =fTrSec[0].GetX(i1);
 412   Double_t xx2=fTrSec[0].GetX(i2);
 413
 414   for (Int_t ns=0; ns<max_sec; ns++) {
 415
 416     Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2];
 417     Int_t nm=fTrSec[ns][i2];
 418     Int_t nu=fTrSec[(ns+1)%max_sec][i2];
 419
 420     AliTRDtimeBin& r1=fTrSec[ns][i1];
 421
 422     for (Int_t is=0; is < r1; is++) {
 423       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
 424
 425       for (Int_t js=0; js < nl+nm+nu; js++) {
 426         const AliTRDcluster *cl;
 427         Double_t x2,   y2,   z2;
 428         Double_t x3=0.,y3=0.;
 429
 430         if (js<nl) {
 431           AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2];
 432           cl=r2[js];
 433           y2=cl->GetY(); z2=cl->GetZ();
 434
 435           x2= xx2*cs+y2*sn;
 436           y2=-xx2*sn+y2*cs;
 437
 438         } else
 439           if (js<nl+nm) {
 440             AliTRDtimeBin& r2=fTrSec[ns][i2];
 441             cl=r2[js-nl];
 442             x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
 443           } else {
 444             AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
 445             cl=r2[js-nl-nm];
 446             y2=cl->GetY(); z2=cl->GetZ();
 447
 448             x2=xx2*cs-y2*sn;
 449             y2=xx2*sn+y2*cs;
 450
 451           }
 452
 453         Double_t zz=z1 - z1/x1*(x1-x2);
 454
 455         if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
 456
 457         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
 458         if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
 459
 460         x[0]=y1;
 461         x[1]=z1;
 462         x[2]=f1trd(x1,y1,x2,y2,x3,y3);
 463
 464         if (TMath::Abs(x[2]) >= fMaxSeedC) continue;
 465
 466         x[3]=f2trd(x1,y1,x2,y2,x3,y3);
 467
 468         if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
 469
 470         x[4]=f3trd(x1,y1,x2,y2,z1,z2);
 471
 472         if (TMath::Abs(x[4]) > fMaxSeedTan) continue;
 473
 474         Double_t a=asin(x[3]);
 475         Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
 476         if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
 477
 478         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
 479         Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2();
 480         Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
 481
 482         Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
 483         Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
 484         Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
 485         Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
 486         Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
 487         Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
 488         Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
 489         Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
 490         Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
 491         Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
 492
 493         c[0]=sy1;
 494         c[1]=0.;       c[2]=sz1;
 495         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
 496         c[6]=f30*sy1;  c[7]=0.;       c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
 497                                       c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
 498         c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
 499         c[13]=f40*sy1*f30+f42*sy2*f32;
 500         c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
 501
 502         UInt_t index=r1.GetIndex(is);
 503         AliTRDtrack *track=new AliTRDtrack(index, x, c, x1, ns*alpha+shift);
 504
 505         Int_t rc=FindProlongation(*track,fTrSec,ns,i2);
 506
 507         if ((rc < 0) ||
 508             (track->GetNclusters() < (i1-i2)*fMinClustersInSeed)) delete track;
 509         else {
 510           fSeeds->AddLast(track); fNseeds++;
 511           cerr<<"found seed "<<fNseeds<<endl;
 512         }
 513       }
 514     }
 515   }
 516
 517   fSeeds->Sort();
 518 }
 519
 520 //___________________________________________________________________
 521 void AliTRDtracker::FindTracks()
 522 {
 523   if (!fClusters) return;
 524
 525   AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
 526   SetUpSectors(fTrSec);
 527
 528   // find tracks
 529
 530   Int_t num_of_time_bins = fTrSec[0].GetNtimeBins();
 531   Int_t nseed=fSeeds->GetEntriesFast();
 532
 533   Int_t nSeedClusters;
 534   for (Int_t i=0; i<nseed; i++) {
 535     cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
 536
 537     AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i));
 538
 539     nSeedClusters = t.GetNclusters();
 540     Double_t alpha=t.GetAlpha();
 541
 542     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 543     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 544     Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
 545
 546     if (FindProlongation(t,fTrSec,ns)) {
 547       cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
 548       if (t.GetNclusters() >= Int_t(fMinClustersInTrack*num_of_time_bins)) {
 549         Int_t label = GetTrackLabel(t);
 550         t.SetLabel(label);
 551         UseClusters(t);
 552
 553         AliTRDtrack *pt = new AliTRDtrack(t);
 554         fTracks->AddLast(pt); fNtracks++;
 555
 556         cerr<<"found track "<<fNtracks<<endl;
 557       }
 558     }
 559     delete fSeeds->RemoveAt(i);
 560   }
 561 }
 562
 563 //__________________________________________________________________
 564 void AliTRDtracker::UseClusters(AliTRDtrack t) {
 565   Int_t ncl=t.GetNclusters();
 566   for (Int_t i=0; i<ncl; i++) {
 567     Int_t index = t.GetClusterIndex(i);
 568     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
 569     c->Use();
 570   }
 571 }
 572
 573 //__________________________________________________________________
 574 Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) {
 575
 576   Int_t label=123456789, index, i, j;
 577   Int_t ncl=t.GetNclusters();
 578   const Int_t range = AliTRDgeometry::kNplan * fGeom->GetTimeMax();
 579   Bool_t label_added;
 580
 581   //  Int_t s[range][2];
 582   Int_t **s = new Int_t* [range];
 583   for (i=0; i<range; i++) {
 584     s[i] = new Int_t[2];
 585   }
 586   for (i=0; i<range; i++) {
 587     s[i][0]=-1;
 588     s[i][1]=0;
 589   }
 590
 591   Int_t t0,t1,t2;
 592   for (i=0; i<ncl; i++) {
 593     index=t.GetClusterIndex(i);
 594     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
 595     t0=c->GetTrackIndex(0);
 596     t1=c->GetTrackIndex(1);
 597     t2=c->GetTrackIndex(2);
 598   }
 599
 600   for (i=0; i<ncl; i++) {
 601     index=t.GetClusterIndex(i);
 602     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
 603     for (Int_t k=0; k<3; k++) {
 604       label=c->GetTrackIndex(k);
 605       label_added=kFALSE; j=0;
 606       if (label >= 0) {
 607         while ( (!label_added) && ( j < range ) ) {
 608           if (s[j][0]==label || s[j][1]==0) {
 609             s[j][0]=label;
 610             s[j][1]=s[j][1]+1;
 611             label_added=kTRUE;
 612           }
 613           j++;
 614         }
 615       }
 616     }
 617   }
 618
 619   Int_t max=0;
 620   label = -123456789;
 621
 622   for (i=0; i<range; i++) {
 623     if (s[i][1]>max) {
 624       max=s[i][1]; label=s[i][0];
 625     }
 626   }
 627   delete []s;
 628   if(max > ncl*fLabelFraction) return label;
 629   else return -1;
 630 }
 631
 632 //___________________________________________________________________
 633
 634 Int_t AliTRDtracker::WriteTracks(const Char_t *filename) {
 635
 636   TDirectory *savedir=gDirectory;
 637
 638   TFile *out=TFile::Open(filename,"RECREATE");
 639
 640   TTree tracktree("TreeT","Tree with TRD tracks");
 641
 642   AliTRDtrack *iotrack=0;
 643   tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
 644
 645   Int_t ntracks=fTracks->GetEntriesFast();
 646
 647   for (Int_t i=0; i<ntracks; i++) {
 648     AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
 649     iotrack=pt;
 650     tracktree.Fill();
 651     cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
 652   }
 653
 654   tracktree.Write();
 655   out->Close();
 656
 657   savedir->cd();
 658
 659   cerr<<"WriteTracks: done"<<endl;
 660   return 0;
 661 }
 662
 663 //_____________________________________________________________________________
 664 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename,
 665 Int_t option)
 666 {
 667   //
 668   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
 669   // from the file. The names of the cluster tree and branches
 670   // should match the ones used in AliTRDclusterizer::WriteClusters()
 671   //
 672
 673   TDirectory *savedir=gDirectory;
 674
 675   TFile *file = TFile::Open(filename);
 676   if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
 677
 678   TTree *tree = (TTree*)file->Get("ClusterTree");
 679   Int_t nentr = (Int_t) tree->GetEntries();
 680   printf("found %d entries in %s.\n",nentr,tree->GetName());
 681
 682   TBranch *branch;
 683   TObjArray *ioArray = new TObjArray(400);
 684
 685   if( option < 0 ) {
 686     //branch = tree->GetBranch("RecPoints");
 687     // changed CBL
 688     branch = tree->GetBranch("TRDrecPoints");
 689
 690     for (Int_t i=0; i<nentr; i++) {
 691       branch->SetAddress(&ioArray);
 692       tree->GetEvent(i);
 693       Int_t npoints = ioArray->GetEntriesFast();
 694       printf("Read %d rec. points from entry %d \n", npoints, i);
 695
 696       for(Int_t j=0; j<npoints; j++) {
 697         AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j);
 698         array->AddLast(p);
 699         ioArray->RemoveAt(j);
 700       }
 701     }
 702   }
 703   else {
 704     branch = tree->GetBranch("TRDcluster");
 705
 706     for (Int_t i=0; i<nentr; i++) {
 707       branch->SetAddress(&ioArray);
 708       tree->GetEvent(i);
 709       Int_t npoints = ioArray->GetEntriesFast();
 710       printf("Read %d clusters from entry %d \n", npoints, i);
 711
 712       for(Int_t j=0; j<npoints; j++) {
 713         AliTRDcluster *c=(AliTRDcluster*)ioArray->UncheckedAt(j);
 714         array->AddLast(c);
 715         ioArray->RemoveAt(j);
 716       }
 717     }
 718   }
 719
 720   file->Close();
 721   savedir->cd();
 722
 723 }
 724