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