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