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