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.14  2001/11/14 10:50:46  cblume
  19 Changes in digits IO. Add merging of summable digits
  20
  21 Revision 1.13  2001/05/30 12:17:47  hristov
  22 Loop variables declared once
  23
  24 Revision 1.12  2001/05/28 17:07:58  hristov
  25 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
  26
  27 Revision 1.8  2000/12/20 13:00:44  cblume
  28 Modifications for the HP-compiler
  29
  30 Revision 1.7  2000/12/08 16:07:02  cblume
  31 Update of the tracking by Sergei
  32
  33 Revision 1.6  2000/11/30 17:38:08  cblume
  34 Changes to get in line with new STEER and EVGEN
  35
  36 Revision 1.5  2000/11/14 14:40:27  cblume
  37 Correction for the Sun compiler (kTRUE and kFALSE)
  38
  39 Revision 1.4  2000/11/10 14:57:52  cblume
  40 Changes in the geometry constants for the DEC compiler
  41
  42 Revision 1.3  2000/10/15 23:40:01  cblume
  43 Remove AliTRDconst
  44
  45 Revision 1.2  2000/10/06 16:49:46  cblume
  46 Made Getters const
  47
  48 Revision 1.1.2.2  2000/10/04 16:34:58  cblume
  49 Replace include files by forward declarations
  50
  51 Revision 1.1.2.1  2000/09/22 14:47:52  cblume
  52 Add the tracking code
  53
  54 */
  55
  56 ////////////////////////////////////////////////////////////////////////////////
  57 //                                                                            //
  58 //  The TRD tracker                                                           //
  59 //                                                                            //
  60 ////////////////////////////////////////////////////////////////////////////////
  61
  62 #include <iostream.h>
  63
  64 #include <TFile.h>
  65 #include <TROOT.h>
  66 #include <TBranch.h>
  67 #include <TTree.h>
  68
  69 #include "AliRun.h"
  70 #include "AliTRD.h"
  71 #include "AliTRDgeometry.h"
  72 #include "AliTRDrecPoint.h"
  73 #include "AliTRDcluster.h"
  74 #include "AliTRDtrack.h"
  75 #include "AliTRDtrackingSector.h"
  76 #include "AliTRDtimeBin.h"
  77
  78 #include "AliTRDtracker.h"
  79
  80 ClassImp(AliTRDtracker)
  81
  82   const  Float_t     AliTRDtracker::fgkSeedDepth          = 0.5;
  83   const  Float_t     AliTRDtracker::fgkSeedStep           = 0.05;
  84   const  Float_t     AliTRDtracker::fgkSeedGap            = 0.25;
  85
  86   const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ12    = 40.;
  87   const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ      = 25.;
  88   const  Float_t     AliTRDtracker::fgkMaxSeedC           = 0.0052;
  89   const  Float_t     AliTRDtracker::fgkMaxSeedTan         = 1.2;
  90   const  Float_t     AliTRDtracker::fgkMaxSeedVertexZ     = 150.;
  91
  92   const  Double_t    AliTRDtracker::fgkSeedErrorSY        = 0.2;
  93   const  Double_t    AliTRDtracker::fgkSeedErrorSY3       = 2.5;
  94   const  Double_t    AliTRDtracker::fgkSeedErrorSZ        = 0.1;
  95
  96   const  Float_t     AliTRDtracker::fgkMinClustersInSeed  = 0.7;
  97
  98   const  Float_t     AliTRDtracker::fgkMinClustersInTrack = 0.5;
  99   const  Float_t     AliTRDtracker::fgkSkipDepth          = 0.05;
 100   const  Float_t     AliTRDtracker::fgkLabelFraction      = 0.5;
 101   const  Float_t     AliTRDtracker::fgkWideRoad           = 20.;
 102
 103   const  Double_t    AliTRDtracker::fgkMaxChi2            = 24.;
 104
 105 //____________________________________________________________________
 106 AliTRDtracker::AliTRDtracker()
 107 {
 108   //
 109   // Default constructor
 110   //
 111
 112   fEvent     = 0;
 113   fGeom      = NULL;
 114
 115   fNclusters = 0;
 116   fClusters  = NULL;
 117   fNseeds    = 0;
 118   fSeeds     = NULL;
 119   fNtracks   = 0;
 120   fTracks    = NULL;
 121
 122   fSY2corr = 0.025;
 123 }
 124
 125 //____________________________________________________________________
 126 AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
 127                   :TNamed(name, title)
 128 {
 129   //
 130   // TRD tracker contructor
 131   //
 132
 133   fEvent     = 0;
 134   fGeom      = NULL;
 135
 136   fNclusters = 0;
 137   fClusters  = new TObjArray(2000);
 138   fNseeds    = 0;
 139   fSeeds     = new TObjArray(20000);
 140   fNtracks   = 0;
 141   fTracks    = new TObjArray(10000);
 142
 143   fSY2corr = 0.025;
 144 }
 145
 146 //___________________________________________________________________
 147 AliTRDtracker::~AliTRDtracker()
 148 {
 149   //
 150   // Destructor
 151   //
 152
 153   delete fClusters;
 154   delete fTracks;
 155   delete fSeeds;
 156   delete fGeom;
 157
 158 }
 159
 160 //___________________________________________________________________
 161 void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd)
 162 {
 163   //
 164   // Do the trackfinding
 165   //
 166
 167   Int_t i;
 168
 169   Int_t inner, outer;
 170   Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
 171   Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
 172   Int_t gap = (Int_t) (fTotalNofTimeBins * fgkSeedGap);
 173   Int_t step = (Int_t) (fTotalNofTimeBins * fgkSeedStep);
 174
 175
 176   //  nSteps = 1;
 177
 178   if (!fClusters) return;
 179
 180   AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
 181   SetUpSectors(fTrSec);
 182
 183   // make a first turn looking for seed ends in the same (n,n)
 184   // and in the adjacent sectors (n,n+1)
 185
 186   for(i=0; i<nSteps; i++) {
 187     printf("step %d out of %d \n", i+1, nSteps);
 188     outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
 189     MakeSeeds(inner,outer, fTrSec, 1, hs, hd);
 190     FindTracks(fTrSec, hs, hd);
 191   }
 192
 193   // make a second turn looking for seed ends in next-to-adjacent
 194   // sectors (n,n+2)
 195
 196   for(i=0; i<nSteps; i++) {
 197     printf("step %d out of %d \n", i+1, nSteps);
 198     outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
 199     MakeSeeds(inner, outer, fTrSec, 2, hs, hd);
 200     FindTracks(fTrSec,hs,hd);
 201   }
 202
 203 }
 204
 205 //_____________________________________________________________________
 206 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const
 207 {
 208   //
 209   // Parametrised "expected" error of the cluster reconstruction in Y
 210   //
 211
 212   Double_t s = 0.08 * 0.08;
 213   return s;
 214
 215 }
 216
 217 //_____________________________________________________________________
 218 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) const
 219 {
 220   //
 221   // Parametrised "expected" error of the cluster reconstruction in Z
 222   //
 223
 224   Double_t s = 6 * 6 /12.;
 225   return s;
 226
 227 }
 228
 229 //_____________________________________________________________________
 230 Double_t f1trd(Double_t x1,Double_t y1,
 231                Double_t x2,Double_t y2,
 232                Double_t x3,Double_t y3)
 233 {
 234   //
 235   // Initial approximation of the track curvature
 236   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 237   //
 238
 239   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 240   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 241                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 242   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 243                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 244
 245   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 246
 247   return -xr*yr/sqrt(xr*xr+yr*yr);
 248
 249 }
 250
 251 //_____________________________________________________________________
 252 Double_t f2trd(Double_t x1,Double_t y1,
 253                Double_t x2,Double_t y2,
 254                Double_t x3,Double_t y3)
 255 {
 256   //
 257   // Initial approximation of the track curvature times center of curvature
 258   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 259   //
 260
 261   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 262   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 263                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 264   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 265                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 266
 267   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 268
 269   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 270
 271 }
 272
 273 //_____________________________________________________________________
 274 Double_t f3trd(Double_t x1,Double_t y1,
 275                Double_t x2,Double_t y2,
 276                Double_t z1,Double_t z2)
 277 {
 278   //
 279   // Initial approximation of the tangent of the track dip angle
 280   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 281   //
 282
 283   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 284
 285 }
 286
 287
 288 //___________________________________________________________________
 289
 290 Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
 291                             Int_t s, Int_t rf, Int_t matchedIndex,
 292                                       TH1F *hs, TH1F *hd)
 293 {
 294   //
 295   // Starting from current position on track=t this function tries
 296   // to extrapolate the track up to timeBin=rf and to confirm prolongation
 297   // if a close cluster is found. *sec is a pointer to allocated
 298   // array of sectors, in which the initial sector has index=s.
 299   //
 300
 301   //  TH1F *hsame = hs;
 302   //  TH1F *hdiff = hd;
 303
 304   //  Bool_t good_match;
 305
 306   const Int_t kTimeBinsToSkip=Int_t(fgkSkipDepth*sec->GetNtimeBins());
 307   Int_t tryAgain=kTimeBinsToSkip;
 308
 309   Double_t alpha=AliTRDgeometry::GetAlpha();
 310
 311   Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
 312
 313   Double_t x0, rho;
 314
 315   for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
 316
 317     Double_t x=sec->GetX(nr);
 318     Double_t ymax=x*TMath::Tan(0.5*alpha);
 319
 320     rho = 0.00295; x0 = 11.09;  // TEC
 321     if(sec->TECframe(nr,t.GetY(),t.GetZ())) {
 322       rho = 1.7; x0 = 33.0;     // G10 frame
 323     }
 324     if(TMath::Abs(x - t.GetX()) > 3) {
 325       rho = 0.0559; x0 = 55.6;  // radiator
 326     }
 327     if (!t.PropagateTo(x,x0,rho,0.139)) {
 328       cerr<<"Can't propagate to x = "<<x<<endl;
 329       return 0;
 330     }
 331
 332     AliTRDtimeBin& timeBin=sec[s][nr];
 333     Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
 334     Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
 335     Double_t road=25.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 336
 337     if (road>fgkWideRoad) {
 338       if (t.GetNclusters() > 4) {
 339         cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
 340       }
 341       return 0;
 342     }
 343
 344     AliTRDcluster *cl=0;
 345     UInt_t index=0;
 346     //    Int_t ncl = 0;
 347
 348     Double_t maxChi2=fgkMaxChi2;
 349
 350     if (timeBin) {
 351
 352       for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
 353         AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
 354
 355         //      good_match = kFALSE;
 356         //      if((c->GetTrackIndex(0) == matchedIndex) ||
 357         //   (c->GetTrackIndex(1) == matchedIndex) ||
 358         //   (c->GetTrackIndex(2) == matchedIndex)) good_match = kTRUE;
 359         //        if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road);
 360         //        if(hdiff) hdiff->Fill(road);
 361
 362         if (c->GetY() > y+road) break;
 363         if (c->IsUsed() > 0) continue;
 364
 365         //      if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z));
 366         //      else hdiff->Fill(TMath::Abs(c->GetZ()-z));
 367
 368         //      if(!good_match) continue;
 369
 370         if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue;
 371
 372         Double_t chi2=t.GetPredictedChi2(c);
 373
 374         //      if((c->GetTrackIndex(0) == matchedIndex) ||
 375         //         (c->GetTrackIndex(1) == matchedIndex) ||
 376         //         (c->GetTrackIndex(2) == matchedIndex))
 377         //        hdiff->Fill(chi2);
 378
 379         //      ncl++;
 380
 381         if (chi2 > maxChi2) continue;
 382         maxChi2=chi2;
 383         cl=c;
 384         index=timeBin.GetIndex(i);
 385       }
 386
 387       if(!cl) {
 388
 389         for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
 390           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
 391
 392           if (c->GetY() > y+road) break;
 393           if (c->IsUsed() > 0) continue;
 394           if((c->GetZ()-z)*(c->GetZ()-z) > 3.25 * 12 * sz2) continue;
 395
 396           Double_t chi2=t.GetPredictedChi2(c);
 397
 398           //      ncl++;
 399
 400           if (chi2 > maxChi2) continue;
 401           maxChi2=chi2;
 402           cl=c;
 403           index=timeBin.GetIndex(i);
 404         }
 405       }
 406
 407     }
 408
 409     if (cl) {
 410
 411       t.Update(cl,maxChi2,index);
 412
 413       tryAgain=kTimeBinsToSkip;
 414     } else {
 415       if (tryAgain==0) break;
 416       if (y > ymax) {
 417         //      cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
 418          s = (s+1) % ns;
 419          if (!t.Rotate(alpha)) {
 420            cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
 421            return 0;
 422          }
 423       } else if (y <-ymax) {
 424         //      cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
 425          s = (s-1+ns) % ns;
 426          if (!t.Rotate(-alpha)) {
 427            cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
 428            return 0;
 429          }
 430       }
 431       if(!sec->TECframe(nr,y,z)) tryAgain--;
 432     }
 433   }
 434
 435   return 1;
 436 }
 437
 438
 439
 440 //_____________________________________________________________________________
 441 void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
 442 {
 443   //
 444   // Opens a ROOT-file with TRD-clusters and reads the cluster-tree in
 445   //
 446
 447   ReadClusters(fClusters, clusterfile);
 448
 449   // get geometry from the file with hits
 450
 451   TFile *fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(hitfile);
 452   if (!fInputFile) {
 453     printf("AliTRDtracker::Open -- ");
 454     printf("Open the ALIROOT-file %s.\n",hitfile);
 455     fInputFile = new TFile(hitfile,"UPDATE");
 456   }
 457   else {
 458     printf("AliTRDtracker::Open -- ");
 459     printf("%s is already open.\n",hitfile);
 460   }
 461
 462   // Get AliRun object from file or create it if not on file
 463
 464   gAlice = (AliRun*) fInputFile->Get("gAlice");
 465   if (gAlice) {
 466     printf("AliTRDtracker::GetEvent -- ");
 467     printf("AliRun object found on file.\n");
 468   }
 469   else {
 470     printf("AliTRDtracker::GetEvent -- ");
 471     printf("Could not find AliRun object.\n");
 472   }
 473
 474   /*
 475   // Import the Trees for the event nEvent in the file
 476   Int_t nparticles = gAlice->GetEvent(fEvent);
 477   cerr<<"nparticles = "<<nparticles<<endl;
 478
 479   if (nparticles <= 0) {
 480     printf("AliTRDtracker::GetEvent -- ");
 481     printf("No entries in the trees for event %d.\n",fEvent);
 482   }
 483   */
 484
 485   AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD");
 486   fGeom = trd->GetGeometry();
 487
 488 }
 489
 490
 491 //_____________________________________________________________________________
 492 void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
 493 {
 494   //
 495   // Fills clusters into TRD tracking_sectors
 496   // Note that the numbering scheme for the TRD tracking_sectors
 497   // differs from that of TRD sectors
 498   //
 499
 500   for (Int_t i=0; i<AliTRDgeometry::Nsect(); i++) sec[i].SetUp();
 501
 502   //  Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
 503
 504   cerr<<"SetUpSectors: sorting clusters"<<endl;
 505
 506   Int_t ncl=fClusters->GetEntriesFast();
 507   UInt_t index;
 508   while (ncl--) {
 509     printf("\r %d left  ",ncl);
 510     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
 511     Int_t detector=c->GetDetector(), localTimeBin=c->GetLocalTimeBin();
 512     Int_t sector=fGeom->GetSector(detector);
 513
 514     Int_t trackingSector = AliTRDgeometry::kNsect - sector - 1;
 515
 516     Int_t tb=sec[sector].GetTimeBin(detector,localTimeBin);
 517     index=ncl;
 518     sec[trackingSector][tb].InsertCluster(c,index);
 519
 520   }
 521   printf("\r\n");
 522 }
 523
 524
 525 //_____________________________________________________________________________
 526 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer,
 527                               AliTRDtrackingSector* fTrSec, Int_t turn,
 528                               TH1F *hs, TH1F *hd)
 529 {
 530   //
 531   // Creates track seeds using clusters in timeBins=i1,i2
 532   //
 533
 534   Int_t i2 = inner, i1 = outer;
 535   Int_t ti[3], to[3];
 536   Int_t nprim = 85600/2;
 537
 538
 539   TH1F *hsame = hs;
 540   TH1F *hdiff = hd;
 541   Bool_t match = false;
 542   Int_t matchedIndex;
 543
 544   // find seeds
 545
 546   Double_t x[5], c[15];
 547   Int_t maxSec=AliTRDgeometry::kNsect;
 548
 549   Double_t alpha=AliTRDgeometry::GetAlpha();
 550   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
 551   Double_t cs=cos(alpha), sn=sin(alpha);
 552   Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
 553
 554   Double_t x1 =fTrSec[0].GetX(i1);
 555   Double_t xx2=fTrSec[0].GetX(i2);
 556
 557
 558   printf("\n");
 559
 560   if((turn != 1)&&(turn != 2)) {
 561     printf("*** Error in MakeSeeds: unexpected turn = %d  \n", turn);
 562     return;
 563   }
 564
 565
 566   for (Int_t ns=0; ns<maxSec; ns++) {
 567
 568     printf("\n MakeSeeds: sector %d \n", ns);
 569
 570     Int_t nl2=fTrSec[(ns-2+maxSec)%maxSec][i2];
 571     Int_t nl=fTrSec[(ns-1+maxSec)%maxSec][i2];
 572     Int_t nm=fTrSec[ns][i2];
 573     Int_t nu=fTrSec[(ns+1)%maxSec][i2];
 574     Int_t nu2=fTrSec[(ns+2)%maxSec][i2];
 575
 576     AliTRDtimeBin& r1=fTrSec[ns][i1];
 577
 578     for (Int_t is=0; is < r1; is++) {
 579       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
 580       for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii);
 581
 582       for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
 583
 584         const AliTRDcluster *cl;
 585         Double_t x2,   y2,   z2;
 586         Double_t x3=0., y3=0.;
 587
 588         if (js<nl2) {
 589           if(turn != 2) continue;
 590           AliTRDtimeBin& r2=fTrSec[(ns-2+maxSec)%maxSec][i2];
 591           cl=r2[js];
 592           y2=cl->GetY(); z2=cl->GetZ();
 593           for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
 594
 595           x2= xx2*cs2+y2*sn2;
 596           y2=-xx2*sn2+y2*cs2;
 597         }
 598         else if (js<nl2+nl) {
 599           if(turn != 1) continue;
 600           AliTRDtimeBin& r2=fTrSec[(ns-1+maxSec)%maxSec][i2];
 601           cl=r2[js-nl2];
 602           y2=cl->GetY(); z2=cl->GetZ();
 603           for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
 604
 605           x2= xx2*cs+y2*sn;
 606           y2=-xx2*sn+y2*cs;
 607
 608         }
 609         else if (js<nl2+nl+nm) {
 610           if(turn != 1) continue;
 611           AliTRDtimeBin& r2=fTrSec[ns][i2];
 612           cl=r2[js-nl2-nl];
 613           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
 614           for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
 615         }
 616         else if (js<nl2+nl+nm+nu) {
 617           if(turn != 1) continue;
 618           AliTRDtimeBin& r2=fTrSec[(ns+1)%maxSec][i2];
 619           cl=r2[js-nl2-nl-nm];
 620           y2=cl->GetY(); z2=cl->GetZ();
 621           for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
 622
 623           x2=xx2*cs-y2*sn;
 624           y2=xx2*sn+y2*cs;
 625
 626         }
 627         else {
 628           if(turn != 2) continue;
 629           AliTRDtimeBin& r2=fTrSec[(ns+2)%maxSec][i2];
 630           cl=r2[js-nl2-nl-nm-nu];
 631           y2=cl->GetY(); z2=cl->GetZ();
 632           for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
 633
 634           x2=xx2*cs2-y2*sn2;
 635           y2=xx2*sn2+y2*cs2;
 636         }
 637
 638
 639         match = false;
 640         matchedIndex = -1;
 641         for (Int_t ii=0; ii<3; ii++) {
 642           // cerr<<"ti["<<ii<<"] = "<<ti[ii]<<"; to["<<ii<<"] = "<<to[ii]<<endl;
 643           if(ti[ii] < 0) continue;
 644           if(ti[ii] >= nprim) continue;
 645           for (Int_t kk=0; kk<3; kk++) {
 646             if(to[kk] < 0) continue;
 647             if(to[kk] >= nprim) continue;
 648             if(ti[ii] == to[kk]) {
 649               //cerr<<"ti["<<ii<<"] = "<<ti[ii]<<" = "<<to[kk]<<" = to["<<kk<<"]"<<endl;
 650               matchedIndex = ti[ii];
 651               match = true;
 652             }
 653           }
 654         }
 655
 656         if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue;
 657
 658         Double_t zz=z1 - z1/x1*(x1-x2);
 659
 660         if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue;
 661
 662         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
 663         if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
 664
 665         x[0]=y1;
 666         x[1]=z1;
 667         x[2]=f1trd(x1,y1,x2,y2,x3,y3);
 668
 669         if (TMath::Abs(x[2]) > fgkMaxSeedC) continue;
 670
 671         x[3]=f2trd(x1,y1,x2,y2,x3,y3);
 672
 673         if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
 674
 675         x[4]=f3trd(x1,y1,x2,y2,z1,z2);
 676
 677         if (TMath::Abs(x[4]) > fgkMaxSeedTan) continue;
 678
 679         Double_t a=asin(x[3]);
 680         Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
 681
 682         if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue;
 683
 684         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
 685         Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2();
 686         Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ;
 687
 688         Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
 689         Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
 690         Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
 691         Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
 692         Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
 693         Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
 694         Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
 695         Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
 696         Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
 697         Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
 698
 699         c[0]=sy1;
 700         c[1]=0.;       c[2]=sz1;
 701         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
 702         c[6]=f30*sy1;  c[7]=0.;       c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
 703                                       c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
 704         c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
 705         c[13]=f40*sy1*f30+f42*sy2*f32;
 706         c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
 707
 708         UInt_t index=r1.GetIndex(is);
 709
 710         AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
 711
 712         Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matchedIndex,hsame,hdiff);
 713
 714         //      if (match) hsame->Fill((Float_t) track->GetNclusters());
 715         //      else hdiff->Fill((Float_t) track->GetNclusters());
 716         //      delete track;
 717         //      continue;
 718
 719         if ((rc < 0) ||
 720             (track->GetNclusters() < (i1-i2)*fgkMinClustersInSeed)) delete track;
 721         else {
 722           fSeeds->AddLast(track); fNseeds++;
 723           printf("\r found seed %d  ", fNseeds);
 724         }
 725       }
 726     }
 727   }
 728
 729   fSeeds->Sort();
 730 }
 731
 732 //_____________________________________________________________________________
 733 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
 734 {
 735   //
 736   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
 737   // from the file. The names of the cluster tree and branches
 738   // should match the ones used in AliTRDclusterizer::WriteClusters()
 739   //
 740
 741   TDirectory *savedir=gDirectory;
 742
 743   TFile *file = TFile::Open(filename);
 744   if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
 745
 746   Char_t treeName[12];
 747   sprintf(treeName,"TreeR%d_TRD",fEvent);
 748   TTree *clusterTree = (TTree*) file->Get(treeName);
 749
 750   TObjArray *clusterArray = new TObjArray(400);
 751
 752   clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
 753
 754   Int_t nEntries = (Int_t) clusterTree->GetEntries();
 755   printf("found %d entries in %s.\n",nEntries,clusterTree->GetName());
 756
 757   // Loop through all entries in the tree
 758   Int_t nbytes;
 759   AliTRDcluster *c = 0;
 760
 761   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
 762
 763     // Import the tree
 764     nbytes += clusterTree->GetEvent(iEntry);
 765
 766     // Get the number of points in the detector
 767     Int_t nCluster = clusterArray->GetEntriesFast();
 768     printf("Read %d clusters from entry %d \n", nCluster, iEntry);
 769
 770     // Loop through all TRD digits
 771     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
 772       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
 773       AliTRDcluster *co = new AliTRDcluster(*c);
 774       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
 775       array->AddLast(co);
 776       delete clusterArray->RemoveAt(iCluster);
 777     }
 778   }
 779
 780   file->Close();
 781   delete clusterArray;
 782   savedir->cd();
 783
 784 }
 785
 786 //___________________________________________________________________
 787 void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd)
 788 {
 789   //
 790   // Finds tracks in TRD
 791   //
 792
 793   TH1F *hsame = hs;
 794   TH1F *hdiff = hd;
 795
 796   Int_t numOfTimeBins = fTrSec[0].GetNtimeBins();
 797   Int_t nseed=fSeeds->GetEntriesFast();
 798
 799   Int_t nSeedClusters;
 800   for (Int_t i=0; i<nseed; i++) {
 801     cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
 802
 803     AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i));
 804
 805     nSeedClusters = t.GetNclusters();
 806     Double_t alpha=t.GetAlpha();
 807
 808     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 809     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 810     Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
 811
 812     Int_t label = GetTrackLabel(t);
 813
 814     if (FindProlongation(t,fTrSec,ns,0,label,hsame,hdiff)) {
 815       cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
 816       if (t.GetNclusters() >= Int_t(fgkMinClustersInTrack*numOfTimeBins)) {
 817         Int_t label = GetTrackLabel(t);
 818         t.SetLabel(label);
 819         t.CookdEdx();
 820         UseClusters(t);
 821
 822         AliTRDtrack *pt = new AliTRDtrack(t);
 823         fTracks->AddLast(pt); fNtracks++;
 824
 825         cerr<<"found track "<<fNtracks<<endl;
 826       }
 827     }
 828     delete fSeeds->RemoveAt(i);
 829     fNseeds--;
 830   }
 831 }
 832
 833 //__________________________________________________________________
 834 void AliTRDtracker::UseClusters(AliTRDtrack t)
 835 {
 836   //
 837   // Mark used cluster
 838   //
 839
 840   Int_t ncl=t.GetNclusters();
 841   for (Int_t i=0; i<ncl; i++) {
 842     Int_t index = t.GetClusterIndex(i);
 843     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
 844     c->Use();
 845   }
 846
 847 }
 848
 849 //__________________________________________________________________
 850 Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t)
 851 {
 852   //
 853   // Get MC label
 854   //
 855
 856   Int_t label=123456789, index, i, j;
 857   Int_t ncl=t.GetNclusters();
 858   const Int_t kRange = AliTRDgeometry::kNplan * fGeom->GetTimeMax();
 859   Bool_t labelAdded;
 860
 861   //  Int_t s[kRange][2];
 862   Int_t **s = new Int_t* [kRange];
 863   for (i=0; i<kRange; i++) {
 864     s[i] = new Int_t[2];
 865   }
 866   for (i=0; i<kRange; i++) {
 867     s[i][0]=-1;
 868     s[i][1]=0;
 869   }
 870
 871   Int_t t0,t1,t2;
 872   for (i=0; i<ncl; i++) {
 873     index=t.GetClusterIndex(i);
 874     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
 875     t0=c->GetTrackIndex(0);
 876     t1=c->GetTrackIndex(1);
 877     t2=c->GetTrackIndex(2);
 878   }
 879
 880   for (i=0; i<ncl; i++) {
 881     index=t.GetClusterIndex(i);
 882     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
 883     for (Int_t k=0; k<3; k++) {
 884       label=c->GetTrackIndex(k);
 885       labelAdded=kFALSE; j=0;
 886       if (label >= 0) {
 887         while ( (!labelAdded) && ( j < kRange ) ) {
 888           if (s[j][0]==label || s[j][1]==0) {
 889             s[j][0]=label;
 890             s[j][1]=s[j][1]+1;
 891             labelAdded=kTRUE;
 892           }
 893           j++;
 894         }
 895       }
 896     }
 897   }
 898
 899   Int_t max=0;
 900   label = -123456789;
 901
 902   for (i=0; i<kRange; i++) {
 903     if (s[i][1]>max) {
 904       max=s[i][1]; label=s[i][0];
 905     }
 906   }
 907   delete []s;
 908   if(max > ncl*fgkLabelFraction) return label;
 909   else return -1;
 910 }
 911
 912 //___________________________________________________________________
 913 Int_t AliTRDtracker::WriteTracks(const Char_t *filename)
 914 {
 915   //
 916   // Write the tracks to the output file
 917   //
 918
 919   TDirectory *savedir=gDirectory;
 920
 921   //TFile *out=TFile::Open(filename,"RECREATE");
 922   TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename);
 923   if (!out) {
 924     printf("AliTRDtracker::Open -- ");
 925     printf("Open the ALIROOT-file %s.\n",filename);
 926     out = new TFile(filename,"RECREATE");
 927   }
 928   else {
 929     printf("AliTRDtracker::Open -- ");
 930     printf("%s is already open.\n",filename);
 931   }
 932
 933   Char_t treeName[12];
 934   sprintf(treeName,"TreeT%d_TRD",fEvent);
 935   TTree tracktree(treeName,"Tree with TRD tracks");
 936
 937   AliTRDtrack *iotrack=0;
 938   tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
 939
 940   Int_t ntracks=fTracks->GetEntriesFast();
 941
 942   for (Int_t i=0; i<ntracks; i++) {
 943     AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
 944     iotrack=pt;
 945     tracktree.Fill();
 946     cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
 947   }
 948
 949   tracktree.Write();
 950   out->Close();
 951
 952   savedir->cd();
 953
 954   cerr<<"WriteTracks: done"<<endl;
 955   return 1;
 956
 957 }
 958