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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  The standard TRD tracker                                                 //
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <Riostream.h>
  25 #include <TFile.h>
  26 #include <TBranch.h>
  27 #include <TTree.h>
  28 #include <TObjArray.h>
  29
  30 #include "AliTRDgeometry.h"
  31 #include "AliTRDparameter.h"
  32 #include "AliTRDgeometryDetail.h"
  33 #include "AliTRDgeometryHole.h"
  34 #include "AliTRDcluster.h"
  35 #include "AliTRDtrack.h"
  36 #include "AliTRDPartID.h"
  37 #include "../TPC/AliTPCtrack.h"
  38
  39 #include "AliTRDtracker.h"
  40
  41 ClassImp(AliTRDtracker)
  42
  43   const  Float_t     AliTRDtracker::fgkSeedDepth          = 0.5;
  44   const  Float_t     AliTRDtracker::fgkSeedStep           = 0.10;
  45   const  Float_t     AliTRDtracker::fgkSeedGap            = 0.25;
  46
  47   const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ12    = 40.;
  48   const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ      = 25.;
  49   const  Float_t     AliTRDtracker::fgkMaxSeedC           = 0.0052;
  50   const  Float_t     AliTRDtracker::fgkMaxSeedTan         = 1.2;
  51   const  Float_t     AliTRDtracker::fgkMaxSeedVertexZ     = 150.;
  52
  53   const  Double_t    AliTRDtracker::fgkSeedErrorSY        = 0.2;
  54   const  Double_t    AliTRDtracker::fgkSeedErrorSY3       = 2.5;
  55   const  Double_t    AliTRDtracker::fgkSeedErrorSZ        = 0.1;
  56
  57   const  Float_t     AliTRDtracker::fgkMinClustersInSeed  = 0.7;
  58
  59   const  Float_t     AliTRDtracker::fgkMinClustersInTrack = 0.5;
  60   const  Float_t     AliTRDtracker::fgkMinFractionOfFoundClusters = 0.8;
  61
  62   const  Float_t     AliTRDtracker::fgkSkipDepth          = 0.3;
  63   const  Float_t     AliTRDtracker::fgkLabelFraction      = 0.8;
  64   const  Float_t     AliTRDtracker::fgkWideRoad           = 20.;
  65
  66   const  Double_t    AliTRDtracker::fgkMaxChi2            = 12.;
  67
  68 const Int_t AliTRDtracker::fgkFirstPlane = 5;
  69 const Int_t AliTRDtracker::fgkLastPlane = 17;
  70
  71
  72 //____________________________________________________________________
  73 AliTRDtracker::AliTRDtracker():AliTracker(),
  74                                fGeom(0),
  75                                fPar(0),
  76                                fNclusters(0),
  77                                fClusters(0),
  78                                fNseeds(0),
  79                                fSeeds(0),
  80                                fNtracks(0),
  81                                fTracks(0),
  82                                fSY2corr(0),
  83                                fSZ2corr(0),
  84                                fTimeBinsPerPlane(0),
  85                                fMaxGap(0),
  86                                fVocal(kFALSE),
  87                                fAddTRDseeds(kFALSE),
  88                                fNoTilt(kFALSE)
  89 {
  90   for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
  91   for(Int_t j=0;j<5;j++)
  92     for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
  93 }
  94 //____________________________________________________________________
  95 AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
  96 {
  97   //
  98   //  Main constructor
  99   //
 100
 101   //Float_t fTzero = 0;
 102
 103   fAddTRDseeds = kFALSE;
 104   fGeom = NULL;
 105   fNoTilt = kFALSE;
 106
 107   TDirectory *savedir=gDirectory;
 108   TFile *in=(TFile*)geomfile;
 109   if (!in->IsOpen()) {
 110     printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
 111     printf("    DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
 112   }
 113   else {
 114     in->cd();
 115 //    in->ls();
 116     fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
 117     fPar  = (AliTRDparameter*) in->Get("TRDparameter");
 118 //    fGeom->Dump();
 119   }
 120
 121   if(fGeom) {
 122     //    fTzero = geo->GetT0();
 123     printf("Found geometry version %d on file \n", fGeom->IsVersion());
 124   }
 125   else {
 126     printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
 127     //printf("The DETAIL TRD geometry will be used\n");
 128     //fGeom = new AliTRDgeometryDetail();
 129     fGeom = new AliTRDgeometryHole();
 130     fGeom->SetPHOShole();
 131     fGeom->SetRICHhole();
 132   }
 133
 134   if (!fPar) {
 135     printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n");
 136     printf("The DEFAULT TRD parameter will be used\n");
 137     fPar = new AliTRDparameter();
 138   }
 139   fPar->ReInit();
 140
 141   savedir->cd();
 142
 143
 144   //  fGeom->SetT0(fTzero);
 145
 146   fNclusters = 0;
 147   fClusters  = new TObjArray(2000);
 148   fNseeds    = 0;
 149   fSeeds     = new TObjArray(2000);
 150   fNtracks   = 0;
 151   fTracks    = new TObjArray(1000);
 152
 153   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 154     Int_t trS = CookSectorIndex(geomS);
 155     fTrSec[trS] = new AliTRDtrackingSector(fGeom, geomS, fPar);
 156     for (Int_t icham=0;icham<AliTRDgeometry::kNcham; icham++){
 157       fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS);
 158     }
 159   }
 160
 161   Float_t tiltAngle = TMath::Abs(fPar->GetTiltingAngle());
 162   if(tiltAngle < 0.1) {
 163     fNoTilt = kTRUE;
 164   }
 165
 166   fSY2corr = 0.2;
 167   fSZ2corr = 120.;
 168
 169   if(fNoTilt && (tiltAngle > 0.1)) fSY2corr = fSY2corr + tiltAngle * 0.05;
 170
 171
 172   // calculate max gap on track
 173
 174   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
 175   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
 176
 177   Double_t dx = (Double_t) fPar->GetTimeBinSize();
 178   Int_t tbAmp = fPar->GetTimeBefore();
 179   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
 180   if(kTRUE) maxAmp = 0;  // intentional until we change the parameter class
 181   Int_t tbDrift = fPar->GetTimeMax();
 182   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
 183
 184   tbDrift = TMath::Min(tbDrift,maxDrift);
 185   tbAmp = TMath::Min(tbAmp,maxAmp);
 186
 187   fTimeBinsPerPlane = tbAmp + tbDrift;
 188   fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fgkSkipDepth);
 189
 190   fVocal = kFALSE;
 191
 192
 193   // Barrel Tracks [SR, 03.04.2003]
 194
 195   fBarrelFile = 0;
 196   fBarrelTree = 0;
 197   fBarrelArray = 0;
 198   fBarrelTrack = 0;
 199
 200   savedir->cd();
 201 }
 202
 203 //___________________________________________________________________
 204 AliTRDtracker::~AliTRDtracker()
 205 {
 206   //
 207   // Destructor of AliTRDtracker
 208   //
 209
 210   if (fClusters) {
 211     fClusters->Delete();
 212     delete fClusters;
 213   }
 214   if (fTracks) {
 215     fTracks->Delete();
 216     delete fTracks;
 217   }
 218   if (fSeeds) {
 219     fSeeds->Delete();
 220     delete fSeeds;
 221   }
 222   delete fGeom;
 223   delete fPar;
 224
 225   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 226     delete fTrSec[geomS];
 227   }
 228 }
 229
 230 //_____________________________________________________________________
 231
 232 void AliTRDtracker::SetBarrelTree(const char *mode) {
 233   //
 234   //
 235   //
 236
 237   if (!IsStoringBarrel()) return;
 238
 239   TDirectory *sav = gDirectory;
 240   if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
 241
 242   char buff[40];
 243   sprintf(buff,  "BarrelTRD_%d_%s", GetEventNumber(), mode);
 244
 245   fBarrelFile->cd();
 246   fBarrelTree = new TTree(buff, "Barrel TPC tracks");
 247
 248   Int_t nRefs = fgkLastPlane - fgkFirstPlane + 1;
 249
 250   if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", nRefs);
 251   for(Int_t i=0; i<nRefs; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
 252
 253   fBarrelTree->Branch("tracks", &fBarrelArray);
 254   sav->cd();
 255 }
 256
 257 //_____________________________________________________________________
 258
 259 void AliTRDtracker::StoreBarrelTrack(AliTRDtrack *ps, Int_t refPlane, Int_t isIn) {
 260   //
 261   //
 262   //
 263
 264   if (!IsStoringBarrel()) return;
 265
 266   static Int_t nClusters;
 267   static Int_t nWrong;
 268   static Double_t chi2;
 269   static Int_t index;
 270   static Bool_t wasLast = kTRUE;
 271
 272   Int_t newClusters, newWrong;
 273   Double_t newChi2;
 274
 275   if (wasLast) {
 276
 277     fBarrelArray->Clear();
 278     nClusters = nWrong = 0;
 279     chi2 = 0.0;
 280     index = 0;
 281     wasLast = kFALSE;
 282   }
 283
 284   fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
 285   ps->GetBarrelTrack(fBarrelTrack);
 286
 287   newClusters = ps->GetNumberOfClusters() - nClusters;
 288   newWrong = ps->GetNWrong() - nWrong;
 289   newChi2 = ps->GetChi2() - chi2;
 290
 291   nClusters =  ps->GetNumberOfClusters();
 292   nWrong = ps->GetNWrong();
 293   chi2 = ps->GetChi2();
 294
 295   if (refPlane != fgkLastPlane) {
 296     fBarrelTrack->SetNClusters(newClusters, newChi2);
 297     fBarrelTrack->SetNWrongClusters(newWrong);
 298   } else {
 299     wasLast = kTRUE;
 300   }
 301
 302   fBarrelTrack->SetRefPlane(refPlane, isIn);
 303 }
 304
 305 //_____________________________________________________________________
 306
 307 Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
 308   //
 309   // Rotates the track when necessary
 310   //
 311
 312   Double_t alpha = AliTRDgeometry::GetAlpha();
 313   Double_t y = track->GetY();
 314   Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
 315
 316   //Int_t ns = AliTRDgeometry::kNsect;
 317   //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
 318
 319   if (y > ymax) {
 320     //s = (s+1) % ns;
 321     if (!track->Rotate(alpha)) return kFALSE;
 322   } else if (y <-ymax) {
 323     //s = (s-1+ns) % ns;
 324     if (!track->Rotate(-alpha)) return kFALSE;
 325   }
 326
 327   return kTRUE;
 328 }
 329
 330 //_____________________________________________________________________
 331 inline Double_t f1trd(Double_t x1,Double_t y1,
 332                       Double_t x2,Double_t y2,
 333                       Double_t x3,Double_t y3)
 334 {
 335   //
 336   // Initial approximation of the track curvature
 337   //
 338   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 339   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 340                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 341   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 342                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 343
 344   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 345
 346   return -xr*yr/sqrt(xr*xr+yr*yr);
 347 }
 348
 349 //_____________________________________________________________________
 350 inline Double_t f2trd(Double_t x1,Double_t y1,
 351                       Double_t x2,Double_t y2,
 352                       Double_t x3,Double_t y3)
 353 {
 354   //
 355   // Initial approximation of the track curvature times X coordinate
 356   // of the center of curvature
 357   //
 358
 359   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 360   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 361                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 362   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 363                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 364
 365   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 366
 367   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 368 }
 369
 370 //_____________________________________________________________________
 371 inline Double_t f3trd(Double_t x1,Double_t y1,
 372                       Double_t x2,Double_t y2,
 373                       Double_t z1,Double_t z2)
 374 {
 375   //
 376   // Initial approximation of the tangent of the track dip angle
 377   //
 378
 379   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 380 }
 381
 382
 383 AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin){
 384   //
 385   //try to find cluster in the backup list
 386   //
 387   AliTRDcluster * cl =0;
 388   UInt_t *indexes = track->GetBackupIndexes();
 389   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
 390     if (indexes[i]==0) break;
 391     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
 392     if (!cli) break;
 393     if (cli->GetLocalTimeBin()!=timebin) continue;
 394     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 395     if (iplane==plane) {
 396       cl = cli;
 397       break;
 398     }
 399   }
 400   return cl;
 401 }
 402
 403
 404 Int_t  AliTRDtracker::GetLastPlane(AliTRDtrack * track){
 405   //
 406   //return last updated plane
 407   Int_t lastplane=0;
 408   UInt_t *indexes = track->GetBackupIndexes();
 409   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
 410     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
 411     if (!cli) break;
 412     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 413     if (iplane>lastplane) {
 414       lastplane = iplane;
 415     }
 416   }
 417   return lastplane;
 418 }
 419 //___________________________________________________________________
 420 Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out)
 421 {
 422   //
 423   // Finds tracks within the TRD. File <inp> is expected to contain seeds
 424   // at the outer part of the TRD. If <inp> is NULL, the seeds
 425   // are found within the TRD if fAddTRDseeds is TRUE.
 426   // The tracks are propagated to the innermost time bin
 427   // of the TRD and stored in file <out>.
 428   //
 429
 430   LoadEvent();
 431
 432   TDirectory *savedir=gDirectory;
 433
 434   char   tname[100];
 435
 436   if (!out->IsOpen()) {
 437     cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n";
 438     return 1;
 439   }
 440
 441   sprintf(tname,"seedTRDtoTPC_%d",GetEventNumber());
 442   TTree tpcTree(tname,"Tree with seeds from TRD at outer TPC pad row");
 443   AliTPCtrack *iotrack=0;
 444   tpcTree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 445
 446   sprintf(tname,"TreeT%d_TRD",GetEventNumber());
 447   TTree trdTree(tname,"TRD tracks at inner TRD time bin");
 448   AliTRDtrack *iotrackTRD=0;
 449   trdTree.Branch("tracks","AliTRDtrack",&iotrackTRD,32000,0);
 450
 451   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 452   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 453
 454   if (inp) {
 455      TFile *in=(TFile*)inp;
 456      if (!in->IsOpen()) {
 457         cerr<<
 458           "AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n";
 459         cerr<<" ... going for seeds finding inside the TRD\n";
 460      }
 461      else {
 462        in->cd();
 463        sprintf(tname,"TRDb_%d",GetEventNumber());
 464        TTree *seedTree=(TTree*)in->Get(tname);
 465        if (!seedTree) {
 466          cerr<<"AliTRDtracker::Clusters2Tracks(): ";
 467          cerr<<"can't get a tree with track seeds !\n";
 468          return 3;
 469        }
 470        AliTRDtrack *seed=new AliTRDtrack;
 471        seedTree->SetBranchAddress("tracks",&seed);
 472
 473        Int_t n=(Int_t)seedTree->GetEntries();
 474        for (Int_t i=0; i<n; i++) {
 475          seedTree->GetEvent(i);
 476          seed->ResetCovariance();
 477          AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
 478          fSeeds->AddLast(tr);
 479          fNseeds++;
 480        }
 481        delete seed;
 482        delete seedTree;
 483      }
 484   }
 485
 486   out->cd();
 487
 488
 489   // find tracks from loaded seeds
 490
 491   Int_t nseed=fSeeds->GetEntriesFast();
 492   Int_t i, found = 0;
 493   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 494
 495   for (i=0; i<nseed; i++) {
 496     AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
 497     FollowProlongation(t, innerTB);
 498     if (t.GetNumberOfClusters() >= foundMin) {
 499       UseClusters(&t);
 500       CookLabel(pt, 1-fgkLabelFraction);
 501       //      t.CookdEdx();
 502     }
 503     iotrackTRD = pt;
 504     trdTree.Fill();
 505     found++;
 506 //    cout<<found<<'\r';
 507
 508     if(PropagateToTPC(t)) {
 509       AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
 510       iotrack = tpc;
 511       tpcTree.Fill();
 512       delete tpc;
 513     }
 514     delete fSeeds->RemoveAt(i);
 515     fNseeds--;
 516   }
 517
 518   cout<<"Number of loaded seeds: "<<nseed<<endl;
 519   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
 520
 521   // after tracks from loaded seeds are found and the corresponding
 522   // clusters are used, look for additional seeds from TRD
 523
 524   if(fAddTRDseeds) {
 525     // Find tracks for the seeds in the TRD
 526     Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 527
 528     Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
 529     Int_t gap = (Int_t) (timeBins * fgkSeedGap);
 530     Int_t step = (Int_t) (timeBins * fgkSeedStep);
 531
 532     // make a first turn with tight cut on initial curvature
 533     for(Int_t turn = 1; turn <= 2; turn++) {
 534       if(turn == 2) {
 535         nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep));
 536         step = (Int_t) (timeBins * (3*fgkSeedStep));
 537       }
 538       for(Int_t i=0; i<nSteps; i++) {
 539         Int_t outer=timeBins-1-i*step;
 540         Int_t inner=outer-gap;
 541
 542         nseed=fSeeds->GetEntriesFast();
 543
 544         MakeSeeds(inner, outer, turn);
 545
 546         nseed=fSeeds->GetEntriesFast();
 547         printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
 548                turn, i, nseed);
 549
 550         for (Int_t i=0; i<nseed; i++) {
 551           AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
 552           FollowProlongation(t,innerTB);
 553           if (t.GetNumberOfClusters() >= foundMin) {
 554             UseClusters(&t);
 555             CookLabel(pt, 1-fgkLabelFraction);
 556             t.CookdEdx();
 557             found++;
 558 //            cout<<found<<'\r';
 559             iotrackTRD = pt;
 560             trdTree.Fill();
 561             if(PropagateToTPC(t)) {
 562               AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
 563               iotrack = tpc;
 564               tpcTree.Fill();
 565               delete tpc;
 566             }
 567           }
 568           delete fSeeds->RemoveAt(i);
 569           fNseeds--;
 570         }
 571       }
 572     }
 573   }
 574   tpcTree.Write();
 575   trdTree.Write();
 576
 577   cout<<"Total number of found tracks: "<<found<<endl;
 578
 579   UnloadEvent();
 580
 581   savedir->cd();
 582
 583   return 0;
 584 }
 585
 586 //___________________________________________________________________
 587 Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
 588 {
 589   //
 590   // Finds tracks within the TRD. The ESD event is expected to contain seeds
 591   // at the outer part of the TRD. The seeds
 592   // are found within the TRD if fAddTRDseeds is TRUE.
 593   // The tracks are propagated to the innermost time bin
 594   // of the TRD and the ESD event is updated
 595   //
 596
 597   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 598   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 599   Int_t nseed = 0;
 600   Int_t found = 0;
 601   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 602
 603   Int_t n = event->GetNumberOfTracks();
 604   for (Int_t i=0; i<n; i++) {
 605     AliESDtrack* seed=event->GetTrack(i);
 606     ULong_t status=seed->GetStatus();
 607     if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
 608     if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
 609     nseed++;
 610
 611     AliTRDtrack* seed2 = new AliTRDtrack(*seed);
 612     //seed2->ResetCovariance();
 613     AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
 614     AliTRDtrack &t=*pt;
 615     FollowProlongation(t, innerTB);
 616     if (t.GetNumberOfClusters() >= foundMin) {
 617       UseClusters(&t);
 618       CookLabel(pt, 1-fgkLabelFraction);
 619       //      t.CookdEdx();
 620     }
 621     found++;
 622 //    cout<<found<<'\r';
 623
 624     if(PropagateToTPC(t)) {
 625       seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
 626     }
 627     delete seed2;
 628     delete pt;
 629   }
 630
 631   cout<<"Number of loaded seeds: "<<nseed<<endl;
 632   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
 633
 634   // after tracks from loaded seeds are found and the corresponding
 635   // clusters are used, look for additional seeds from TRD
 636
 637   if(fAddTRDseeds) {
 638     // Find tracks for the seeds in the TRD
 639     Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 640
 641     Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
 642     Int_t gap = (Int_t) (timeBins * fgkSeedGap);
 643     Int_t step = (Int_t) (timeBins * fgkSeedStep);
 644
 645     // make a first turn with tight cut on initial curvature
 646     for(Int_t turn = 1; turn <= 2; turn++) {
 647       if(turn == 2) {
 648         nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep));
 649         step = (Int_t) (timeBins * (3*fgkSeedStep));
 650       }
 651       for(Int_t i=0; i<nSteps; i++) {
 652         Int_t outer=timeBins-1-i*step;
 653         Int_t inner=outer-gap;
 654
 655         nseed=fSeeds->GetEntriesFast();
 656
 657         MakeSeeds(inner, outer, turn);
 658
 659         nseed=fSeeds->GetEntriesFast();
 660         printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
 661                turn, i, nseed);
 662
 663         for (Int_t i=0; i<nseed; i++) {
 664           AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
 665           FollowProlongation(t,innerTB);
 666           if (t.GetNumberOfClusters() >= foundMin) {
 667             UseClusters(&t);
 668             CookLabel(pt, 1-fgkLabelFraction);
 669             t.CookdEdx();
 670             found++;
 671 //            cout<<found<<'\r';
 672             if(PropagateToTPC(t)) {
 673               AliESDtrack track;
 674               track.UpdateTrackParams(pt,AliESDtrack::kTRDin);
 675               event->AddTrack(&track);
 676             }
 677           }
 678           delete fSeeds->RemoveAt(i);
 679           fNseeds--;
 680         }
 681       }
 682     }
 683   }
 684
 685   cout<<"Total number of found tracks: "<<found<<endl;
 686
 687   return 0;
 688 }
 689
 690
 691
 692 //_____________________________________________________________________________
 693 Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
 694   //
 695   // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and
 696   // backpropagated by the TPC tracker. Each seed is first propagated
 697   // to the TRD, and then its prolongation is searched in the TRD.
 698   // If sufficiently long continuation of the track is found in the TRD
 699   // the track is updated, otherwise it's stored as originaly defined
 700   // by the TPC tracker.
 701   //
 702
 703   LoadEvent();
 704
 705   TDirectory *savedir=gDirectory;
 706
 707   TFile *in=(TFile*)inp;
 708
 709   if (!in->IsOpen()) {
 710      cerr<<"AliTRDtracker::PropagateBack(): ";
 711      cerr<<"file with back propagated TPC tracks is not open !\n";
 712      return 1;
 713   }
 714
 715   if (!out->IsOpen()) {
 716      cerr<<"AliTRDtracker::PropagateBack(): ";
 717      cerr<<"file for back propagated TRD tracks is not open !\n";
 718      return 2;
 719   }
 720
 721   in->cd();
 722   char   tname[100];
 723   sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber());
 724   TTree *seedTree=(TTree*)in->Get(tname);
 725   if (!seedTree) {
 726      cerr<<"AliTRDtracker::PropagateBack(): ";
 727      cerr<<"can't get a tree with seeds from TPC !\n";
 728      cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n";
 729      return 3;
 730   }
 731
 732   AliTPCtrack *seed=new AliTPCtrack;
 733   seedTree->SetBranchAddress("tracks",&seed);
 734
 735   Int_t n=(Int_t)seedTree->GetEntries();
 736   for (Int_t i=0; i<n; i++) {
 737      seedTree->GetEvent(i);
 738      Int_t lbl = seed->GetLabel();
 739      AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
 740      tr->SetSeedLabel(lbl);
 741      fSeeds->AddLast(tr);
 742      fNseeds++;
 743   }
 744
 745   delete seed;
 746   delete seedTree;
 747
 748   out->cd();
 749
 750   AliTPCtrack *otrack=0;
 751
 752   sprintf(tname,"seedsTRDtoTOF1_%d",GetEventNumber());
 753   TTree tofTree1(tname,"Tracks back propagated through TPC and TRD");
 754   tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 755
 756   sprintf(tname,"seedsTRDtoTOF2_%d",GetEventNumber());
 757   TTree tofTree2(tname,"Tracks back propagated through TPC and TRD");
 758   tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 759
 760   sprintf(tname,"seedsTRDtoPHOS_%d",GetEventNumber());
 761   TTree phosTree(tname,"Tracks back propagated through TPC and TRD");
 762   phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 763
 764   sprintf(tname,"seedsTRDtoRICH_%d",GetEventNumber());
 765   TTree richTree(tname,"Tracks back propagated through TPC and TRD");
 766   richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 767
 768   sprintf(tname,"TRDb_%d",GetEventNumber());
 769   TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin");
 770   AliTRDtrack *otrackTRD=0;
 771   trdTree.Branch("tracks","AliTRDtrack",&otrackTRD,32000,0);
 772
 773   if (IsStoringBarrel()) SetBarrelTree("back");
 774   out->cd();
 775
 776   Int_t found=0;
 777   Int_t nseed=fSeeds->GetEntriesFast();
 778
 779   //  Float_t foundMin = fgkMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan();
 780   Float_t foundMin = 40;
 781
 782   Int_t outermostTB  = fTrSec[0]->GetOuterTimeBin();
 783
 784   for (Int_t i=0; i<nseed; i++) {
 785
 786     AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
 787     Int_t expectedClr = FollowBackProlongation(s);
 788
 789     if (IsStoringBarrel()) {
 790       StoreBarrelTrack(ps, fgkLastPlane, kTrackBack);
 791       fBarrelTree->Fill();
 792     }
 793
 794     Int_t foundClr = s.GetNumberOfClusters();
 795     Int_t lastTB = fTrSec[0]->GetLayerNumber(s.GetX());
 796
 797     //    printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
 798     //     i, foundClr, expectedClr, foundMin);
 799
 800     if (foundClr >= foundMin) {
 801       if(foundClr >= 2) {
 802         s.CookdEdx();
 803         CookLabel(ps, 1-fgkLabelFraction);
 804         UseClusters(ps);
 805       }
 806
 807       // Propagate to outer reference plane [SR, GSI, 18.02.2003]
 808       ps->PropagateTo(364.8);
 809       otrackTRD=ps;
 810       trdTree.Fill();
 811       found++;
 812 //      cout<<found<<'\r';
 813     }
 814
 815     if(((expectedClr < 10) && (lastTB == outermostTB)) ||
 816        ((expectedClr >= 10) &&
 817         (((Float_t) foundClr) / ((Float_t) expectedClr) >=
 818          fgkMinFractionOfFoundClusters) && (lastTB == outermostTB))) {
 819
 820       Double_t xTOF = 375.5;
 821
 822       if(PropagateToOuterPlane(s,xTOF)) {
 823         AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 824         otrack = pt;
 825         tofTree1.Fill();
 826         delete pt;
 827
 828         xTOF = 381.5;
 829
 830         if(PropagateToOuterPlane(s,xTOF)) {
 831           AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 832           otrack = pt;
 833           tofTree2.Fill();
 834           delete pt;
 835
 836           Double_t xPHOS = 460.;
 837
 838           if(PropagateToOuterPlane(s,xPHOS)) {
 839             AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 840             otrack = pt;
 841             phosTree.Fill();
 842             delete pt;
 843
 844             Double_t xRICH = 490+1.267;
 845
 846             if(PropagateToOuterPlane(s,xRICH)) {
 847               AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 848               otrack = pt;
 849               richTree.Fill();
 850               delete pt;
 851             }
 852           }
 853         }
 854       }
 855     }
 856   }
 857
 858
 859   out->cd();
 860   tofTree1.Write();
 861   tofTree2.Write();
 862   phosTree.Write();
 863   richTree.Write();
 864   trdTree.Write();
 865
 866   if (IsStoringBarrel()) { // [SR, 03.04.2003]
 867     fBarrelFile->cd();
 868     fBarrelTree->Write();
 869     fBarrelFile->Flush();
 870   }
 871
 872   savedir->cd();
 873   cerr<<"Number of seeds: "<<nseed<<endl;
 874   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
 875
 876   UnloadEvent();
 877
 878   return 0;
 879
 880 }
 881
 882 //_____________________________________________________________________________
 883 Int_t AliTRDtracker::PropagateBack(AliESD* event) {
 884   //
 885   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
 886   // backpropagated by the TPC tracker. Each seed is first propagated
 887   // to the TRD, and then its prolongation is searched in the TRD.
 888   // If sufficiently long continuation of the track is found in the TRD
 889   // the track is updated, otherwise it's stored as originaly defined
 890   // by the TPC tracker.
 891   //
 892
 893   Int_t found=0;
 894   Float_t foundMin = 40;
 895
 896   Int_t n = event->GetNumberOfTracks();
 897   for (Int_t i=0; i<n; i++) {
 898     AliESDtrack* seed=event->GetTrack(i);
 899     ULong_t status=seed->GetStatus();
 900     if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
 901     if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
 902
 903     Int_t lbl = seed->GetLabel();
 904     AliTRDtrack *track = new AliTRDtrack(*seed);
 905     track->SetSeedLabel(lbl);
 906     fNseeds++;
 907
 908     Int_t expectedClr = FollowBackProlongation(*track);
 909     if (track->GetNumberOfClusters()<expectedClr/3){
 910       AliTRDtrack *track1 = new AliTRDtrack(*seed);
 911       track1->SetSeedLabel(lbl);
 912       FollowBackProlongation(*track1);
 913       AliTRDtrack *track2= new AliTRDtrack(*seed);
 914       track->SetSeedLabel(lbl);
 915       FollowBackProlongation(*track2);
 916       delete track1;
 917       delete track2;
 918     }
 919
 920     Int_t foundClr = track->GetNumberOfClusters();
 921     if (foundClr >= foundMin) {
 922       if(foundClr >= 2) {
 923         track->CookdEdx();
 924 //      CookLabel(track, 1-fgkLabelFraction);
 925         UseClusters(track);
 926       }
 927
 928       // Propagate to outer reference plane [SR, GSI, 18.02.2003]
 929 //      track->PropagateTo(364.8);  why?
 930
 931       //seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 932       //found++;
 933     }
 934
 935     //Propagation to the TOF (I.Belikov)
 936
 937     if (track->GetStop()==kFALSE){
 938       Double_t xTOF = 375.5;
 939       PropagateToOuterPlane(*track,xTOF);
 940
 941       Double_t xtof=378.;
 942       Double_t c2=track->GetC()*xtof - track->GetEta();
 943       if (TMath::Abs(c2)>=0.9999999) continue;
 944
 945       Double_t ymax=xtof*TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
 946       Double_t y=track->GetYat(xtof);
 947       if (y > ymax) {
 948         if (!track->Rotate(AliTRDgeometry::GetAlpha())) return 1;
 949       } else if (y <-ymax) {
 950         if (!track->Rotate(-AliTRDgeometry::GetAlpha())) return 1;
 951       }
 952
 953       if (track->PropagateTo(xtof)) {
 954         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 955         if (track->GetNumberOfClusters()>foundMin) found++;
 956       }
 957     }else{
 958       if (track->GetNumberOfClusters()>15&&track->GetNumberOfClusters()>0.5*expectedClr){
 959         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 960         seed->UpdateTrackParams(track, AliESDtrack::kTRDStop);
 961         found++;
 962       }
 963     }
 964
 965     delete track;
 966
 967     //End of propagation to the TOF
 968     //if (foundClr>foundMin)
 969     //  seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 970
 971
 972   }
 973
 974   cerr<<"Number of seeds: "<<fNseeds<<endl;
 975   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
 976
 977   fSeeds->Clear(); fNseeds=0;
 978
 979   return 0;
 980
 981 }
 982
 983 //_____________________________________________________________________________
 984 Int_t AliTRDtracker::RefitInward(AliESD* event)
 985 {
 986   //
 987   // Refits tracks within the TRD. The ESD event is expected to contain seeds
 988   // at the outer part of the TRD.
 989   // The tracks are propagated to the innermost time bin
 990   // of the TRD and the ESD event is updated
 991   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
 992   //
 993
 994   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 995   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 996   Int_t nseed = 0;
 997   Int_t found = 0;
 998   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 999
1000   Int_t n = event->GetNumberOfTracks();
1001   for (Int_t i=0; i<n; i++) {
1002     AliESDtrack* seed=event->GetTrack(i);
1003     ULong_t status=seed->GetStatus();
1004     if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
1005     if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
1006     nseed++;
1007
1008     AliTRDtrack* seed2 = new AliTRDtrack(*seed);
1009     seed2->ResetCovariance(5.);
1010     AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
1011     UInt_t * indexes2 = seed2->GetIndexes();
1012     UInt_t * indexes3 = pt->GetBackupIndexes();
1013     for (Int_t i=0;i<200;i++) {
1014       if (indexes2[i]==0) break;
1015       indexes3[i] = indexes2[i];
1016     }
1017     //AliTRDtrack *pt = seed2;
1018     AliTRDtrack &t=*pt;
1019     FollowProlongation(t, innerTB);
1020
1021     if (t.GetNumberOfClusters()<seed->GetTRDclusters(indexes3)*0.5){
1022       // debug  - why we dont go back?
1023       AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
1024       UInt_t * indexes2 = seed2->GetIndexes();
1025       UInt_t * indexes3 = pt2->GetBackupIndexes();
1026       for (Int_t i=0;i<200;i++) {
1027         if (indexes2[i]==0) break;
1028         indexes3[i] = indexes2[i];
1029       }
1030       FollowProlongation(*pt2, innerTB);
1031       delete pt2;
1032     }
1033
1034     if (t.GetNumberOfClusters() >= foundMin) {
1035       //      UseClusters(&t);
1036       //CookLabel(pt, 1-fgkLabelFraction);
1037       //      t.CookdEdx();
1038     }
1039     found++;
1040 //    cout<<found<<'\r';
1041
1042     if(PropagateToTPC(t)) {
1043       seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
1044     }
1045     delete seed2;
1046     delete pt;
1047   }
1048
1049   cout<<"Number of loaded seeds: "<<nseed<<endl;
1050   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
1051
1052   return 0;
1053
1054 }
1055
1056
1057 //---------------------------------------------------------------------------
1058 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
1059 {
1060   // Starting from current position on track=t this function tries
1061   // to extrapolate the track up to timeBin=0 and to confirm prolongation
1062   // if a close cluster is found. Returns the number of clusters
1063   // expected to be found in sensitive layers
1064
1065   Float_t  wIndex, wTB, wChi2;
1066   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
1067   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
1068   Float_t  wPx, wPy, wPz, wC;
1069   Double_t px, py, pz;
1070   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
1071   Int_t lastplane = GetLastPlane(&t);
1072
1073   Int_t trackIndex = t.GetLabel();
1074
1075   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1076
1077   Int_t tryAgain=fMaxGap;
1078
1079   Double_t alpha=t.GetAlpha();
1080   alpha = TVector2::Phi_0_2pi(alpha);
1081
1082   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1083   Double_t radLength, rho, x, dx, y, ymax, z;
1084
1085   Int_t expectedNumberOfClusters = 0;
1086   Bool_t lookForCluster;
1087
1088   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1089
1090
1091   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
1092
1093     y = t.GetY(); z = t.GetZ();
1094
1095     // first propagate to the inner surface of the current time bin
1096     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1097     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1098     if(!t.PropagateTo(x,radLength,rho)) break;
1099     y = t.GetY();
1100     ymax = x*TMath::Tan(0.5*alpha);
1101     if (y > ymax) {
1102       s = (s+1) % ns;
1103       if (!t.Rotate(alpha)) break;
1104       if(!t.PropagateTo(x,radLength,rho)) break;
1105     } else if (y <-ymax) {
1106       s = (s-1+ns) % ns;
1107       if (!t.Rotate(-alpha)) break;
1108       if(!t.PropagateTo(x,radLength,rho)) break;
1109     }
1110
1111     y = t.GetY(); z = t.GetZ();
1112
1113     // now propagate to the middle plane of the next time bin
1114     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1115     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1116     if(!t.PropagateTo(x,radLength,rho)) break;
1117     y = t.GetY();
1118     ymax = x*TMath::Tan(0.5*alpha);
1119     if (y > ymax) {
1120       s = (s+1) % ns;
1121       if (!t.Rotate(alpha)) break;
1122       if(!t.PropagateTo(x,radLength,rho)) break;
1123     } else if (y <-ymax) {
1124       s = (s-1+ns) % ns;
1125       if (!t.Rotate(-alpha)) break;
1126       if(!t.PropagateTo(x,radLength,rho)) break;
1127     }
1128
1129
1130     if(lookForCluster) {
1131
1132       expectedNumberOfClusters++;
1133       wIndex = (Float_t) t.GetLabel();
1134       wTB = nr;
1135
1136       AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr-1));
1137
1138       Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
1139       Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
1140
1141       Double_t road;
1142       if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
1143       else return expectedNumberOfClusters;
1144
1145       wYrt = (Float_t) y;
1146       wZrt = (Float_t) z;
1147       wYwindow = (Float_t) road;
1148       t.GetPxPyPz(px,py,pz);
1149       wPx = (Float_t) px;
1150       wPy = (Float_t) py;
1151       wPz = (Float_t) pz;
1152       wC  = (Float_t) t.GetC();
1153       wSigmaC2 = (Float_t) t.GetSigmaC2();
1154       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
1155       wSigmaY2 = (Float_t) t.GetSigmaY2();
1156       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
1157       wChi2 = -1;
1158
1159
1160       AliTRDcluster *cl=0;
1161       UInt_t index=0;
1162
1163       Double_t maxChi2=fgkMaxChi2;
1164
1165       wYclosest = 12345678;
1166       wYcorrect = 12345678;
1167       wZclosest = 12345678;
1168       wZcorrect = 12345678;
1169       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);
1170
1171       // Find the closest correct cluster for debugging purposes
1172       if (timeBin) {
1173         Float_t minDY = 1000000;
1174         for (Int_t i=0; i<timeBin; i++) {
1175           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1176           if((c->GetLabel(0) != trackIndex) &&
1177              (c->GetLabel(1) != trackIndex) &&
1178              (c->GetLabel(2) != trackIndex)) continue;
1179           if(TMath::Abs(c->GetY() - y) > minDY) continue;
1180           minDY = TMath::Abs(c->GetY() - y);
1181           wYcorrect = c->GetY();
1182           wZcorrect = c->GetZ();
1183
1184           Double_t h01 = GetTiltFactor(c);
1185           wChi2 = t.GetPredictedChi2(c, h01);
1186         }
1187       }
1188
1189       // Now go for the real cluster search
1190
1191       if (timeBin) {
1192         //
1193         //find cluster in history
1194         cl =0;
1195
1196         AliTRDcluster * cl0 = timeBin[0];
1197         if (!cl0) {
1198           continue;
1199         }
1200         Int_t plane = fGeom->GetPlane(cl0->GetDetector());
1201         if (plane>lastplane) continue;
1202         Int_t timebin = cl0->GetLocalTimeBin();
1203         AliTRDcluster * cl2= GetCluster(&t,plane, timebin);
1204         if (cl2) {
1205           cl =cl2;
1206           Double_t h01 = GetTiltFactor(cl);
1207           maxChi2=t.GetPredictedChi2(cl,h01);
1208         }
1209         if ((!cl) && road>fgkWideRoad) {
1210           //if (t.GetNumberOfClusters()>4)
1211           //  cerr<<t.GetNumberOfClusters()
1212           //    <<"FindProlongation warning: Too broad road !\n";
1213           continue;
1214         }
1215
1216
1217         if(!cl){
1218
1219           for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1220             AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1221             if (c->GetY() > y+road) break;
1222             if (c->IsUsed() > 0) continue;
1223             if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1224
1225             Double_t h01 = GetTiltFactor(c);
1226             Double_t chi2=t.GetPredictedChi2(c,h01);
1227
1228             if (chi2 > maxChi2) continue;
1229             maxChi2=chi2;
1230             cl=c;
1231             index=timeBin.GetIndex(i);
1232           }
1233         }
1234
1235         if(!cl) {
1236
1237           for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1238             AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1239
1240             if (c->GetY() > y+road) break;
1241             if (c->IsUsed() > 0) continue;
1242             if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue;
1243
1244             Double_t h01 = GetTiltFactor(c);
1245             Double_t chi2=t.GetPredictedChi2(c, h01);
1246
1247             if (chi2 > maxChi2) continue;
1248             maxChi2=chi2;
1249             cl=c;
1250             index=timeBin.GetIndex(i);
1251           }
1252         }
1253         if (cl) {
1254           wYclosest = cl->GetY();
1255           wZclosest = cl->GetZ();
1256           Double_t h01 = GetTiltFactor(cl);
1257
1258           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1259           //printf("Track   position\t%f\t%f\t%f\n",t.GetX(),t.GetY(),t.GetZ());
1260           //printf("Cluster position\t%d\t%f\t%f\n",cl->GetLocalTimeBin(),cl->GetY(),cl->GetZ());
1261           Int_t det = cl->GetDetector();
1262           Int_t plane = fGeom->GetPlane(det);
1263
1264           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1265             //if(!t.Update(cl,maxChi2,index,h01)) {
1266             //if(!tryAgain--) return 0;
1267           }
1268           else tryAgain=fMaxGap;
1269         }
1270         else {
1271           //if (tryAgain==0) break;
1272           tryAgain--;
1273         }
1274
1275         /*
1276         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1277
1278           printf(" %f", wIndex);       //1
1279           printf(" %f", wTB);          //2
1280           printf(" %f", wYrt);         //3
1281           printf(" %f", wYclosest);    //4
1282           printf(" %f", wYcorrect);    //5
1283           printf(" %f", wYwindow);     //6
1284           printf(" %f", wZrt);         //7
1285           printf(" %f", wZclosest);    //8
1286           printf(" %f", wZcorrect);    //9
1287           printf(" %f", wZwindow);     //10
1288           printf(" %f", wPx);          //11
1289           printf(" %f", wPy);          //12
1290           printf(" %f", wPz);          //13
1291           printf(" %f", wSigmaC2*1000000);  //14
1292           printf(" %f", wSigmaTgl2*1000);   //15
1293           printf(" %f", wSigmaY2);     //16
1294           //      printf(" %f", wSigmaZ2);     //17
1295           printf(" %f", wChi2);     //17
1296           printf(" %f", wC);           //18
1297           printf("\n");
1298         }
1299         */
1300       }
1301     }
1302   }
1303   return expectedNumberOfClusters;
1304
1305
1306 }
1307
1308 //___________________________________________________________________
1309
1310 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
1311 {
1312   // Starting from current radial position of track <t> this function
1313   // extrapolates the track up to outer timebin and in the sensitive
1314   // layers confirms prolongation if a close cluster is found.
1315   // Returns the number of clusters expected to be found in sensitive layers
1316
1317   Float_t  wIndex, wTB, wChi2;
1318   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
1319   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
1320   Float_t  wPx, wPy, wPz, wC;
1321   Double_t px, py, pz;
1322   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
1323
1324   Int_t trackIndex = t.GetLabel();
1325   Int_t tryAgain=fMaxGap;
1326
1327   Double_t alpha=t.GetAlpha();
1328   TVector2::Phi_0_2pi(alpha);
1329
1330   Int_t s;
1331
1332   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
1333   Double_t radLength, rho, x, dx, y, ymax = 0, z;
1334   Bool_t lookForCluster;
1335
1336   Int_t expectedNumberOfClusters = 0;
1337   x = t.GetX();
1338
1339   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1340
1341   Int_t nRefPlane = fgkFirstPlane;
1342   Bool_t isNewLayer = kFALSE;
1343
1344   Double_t chi2;
1345   Double_t minDY;
1346   Int_t zone =-10;
1347   Int_t nr;
1348   for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) {
1349
1350     y = t.GetY();
1351     z = t.GetZ();
1352
1353     // first propagate to the outer surface of the current time bin
1354
1355     s = t.GetSector();
1356     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1357     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2;
1358     y = t.GetY();
1359     z = t.GetZ();
1360
1361     if(!t.PropagateTo(x,radLength,rho)) break;
1362     //    if (!AdjustSector(&t)) break;
1363     //
1364     // MI -fix untill correct material desription will be implemented
1365     //
1366     Float_t angle =  t.GetAlpha();  // MI - if rotation - we go through the material -
1367     if (!AdjustSector(&t)) break;
1368     if (TMath::Abs(angle -  t.GetAlpha())>0.000001) break; //better to stop track
1369     Int_t currentzone = fTrSec[s]->GetLayer(nr)->GetZone(z);
1370     if (currentzone==-10) break;  // we are in the frame
1371     if (currentzone>-10){   // layer knows where we are
1372       if (zone==-10) zone = currentzone;
1373       if (zone!=currentzone) break;
1374     }
1375     //
1376     //
1377     s = t.GetSector();
1378     if (!t.PropagateTo(x,radLength,rho)) break;
1379
1380     y = t.GetY();
1381     z = t.GetZ();
1382
1383     // Barrel Tracks [SR, 04.04.2003]
1384
1385     s = t.GetSector();
1386     if (fTrSec[s]->GetLayer(nr)->IsSensitive() !=
1387         fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) {
1388
1389 //      if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
1390     }
1391
1392     if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() &&
1393           ! fTrSec[s]->GetLayer(nr)->IsSensitive()) {
1394       isNewLayer = kTRUE;
1395     } else {isNewLayer = kFALSE;}
1396
1397     y = t.GetY();
1398     z = t.GetZ();
1399
1400     // now propagate to the middle plane of the next time bin
1401     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1402
1403     x = fTrSec[s]->GetLayer(nr+1)->GetX();
1404       if(!t.PropagateTo(x,radLength,rho)) break;
1405     if (!AdjustSector(&t)) break;
1406     s = t.GetSector();
1407       if(!t.PropagateTo(x,radLength,rho)) break;
1408
1409     y = t.GetY();
1410     z = t.GetZ();
1411
1412     if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
1413     //    printf("label %d, pl %d, lookForCluster %d \n",
1414     //     trackIndex, nr+1, lookForCluster);
1415
1416     if(lookForCluster) {
1417       expectedNumberOfClusters++;
1418
1419       wIndex = (Float_t) t.GetLabel();
1420       wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
1421
1422       AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr+1));
1423       Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
1424       Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
1425       if((t.GetSigmaY2() + sy2) < 0) break;
1426       Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2);
1427       Double_t y=t.GetY(), z=t.GetZ();
1428
1429       wYrt = (Float_t) y;
1430       wZrt = (Float_t) z;
1431       wYwindow = (Float_t) road;
1432       t.GetPxPyPz(px,py,pz);
1433       wPx = (Float_t) px;
1434       wPy = (Float_t) py;
1435       wPz = (Float_t) pz;
1436       wC  = (Float_t) t.GetC();
1437       wSigmaC2 = (Float_t) t.GetSigmaC2();
1438       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
1439       wSigmaY2 = (Float_t) t.GetSigmaY2();
1440       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
1441       wChi2 = -1;
1442
1443       if (road>fgkWideRoad) {
1444         if (t.GetNumberOfClusters()>4)
1445           cerr<<t.GetNumberOfClusters()
1446               <<"FindProlongation warning: Too broad road !\n";
1447         return 0;
1448       }
1449
1450       AliTRDcluster *cl=0;
1451       UInt_t index=0;
1452
1453       Double_t maxChi2=fgkMaxChi2;
1454
1455       if (isNewLayer) {
1456         road = 3 * road;
1457         //sz2 = 3 * sz2;
1458         maxChi2 = 10 * fgkMaxChi2;
1459       }
1460
1461       if (nRefPlane == fgkFirstPlane) maxChi2 = 20 * fgkMaxChi2;
1462       if (nRefPlane == fgkFirstPlane+2) maxChi2 = 15 * fgkMaxChi2;
1463       if (t.GetNRotate() > 0) maxChi2 = 3 * maxChi2;
1464
1465
1466       wYclosest = 12345678;
1467       wYcorrect = 12345678;
1468       wZclosest = 12345678;
1469       wZcorrect = 12345678;
1470       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);
1471
1472       // Find the closest correct cluster for debugging purposes
1473       if (timeBin) {
1474         minDY = 1000000;
1475         for (Int_t i=0; i<timeBin; i++) {
1476           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1477           if((c->GetLabel(0) != trackIndex) &&
1478              (c->GetLabel(1) != trackIndex) &&
1479              (c->GetLabel(2) != trackIndex)) continue;
1480           if(TMath::Abs(c->GetY() - y) > minDY) continue;
1481           //minDY = TMath::Abs(c->GetY() - y);
1482           minDY = c->GetY() - y;
1483           wYcorrect = c->GetY();
1484           wZcorrect = c->GetZ();
1485
1486           Double_t h01 = GetTiltFactor(c);
1487           wChi2 = t.GetPredictedChi2(c, h01);
1488         }
1489       }
1490
1491       // Now go for the real cluster search
1492
1493       if (timeBin) {
1494
1495         for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1496           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1497           if (c->GetY() > y+road) break;
1498           if (c->IsUsed() > 0) continue;
1499           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1500
1501           Double_t h01 = GetTiltFactor(c);
1502           chi2=t.GetPredictedChi2(c,h01);
1503
1504           if (chi2 > maxChi2) continue;
1505           maxChi2=chi2;
1506           cl=c;
1507           index=timeBin.GetIndex(i);
1508
1509           //check is correct
1510           if((c->GetLabel(0) != trackIndex) &&
1511              (c->GetLabel(1) != trackIndex) &&
1512              (c->GetLabel(2) != trackIndex)) t.AddNWrong();
1513         }
1514
1515         if(!cl) {
1516
1517           for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1518             AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1519
1520             if (c->GetY() > y+road) break;
1521             if (c->IsUsed() > 0) continue;
1522             if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
1523
1524             Double_t h01 = GetTiltFactor(c);
1525             chi2=t.GetPredictedChi2(c,h01);
1526
1527             if (chi2 > maxChi2) continue;
1528             maxChi2=chi2;
1529             cl=c;
1530             index=timeBin.GetIndex(i);
1531           }
1532         }
1533
1534         if (cl) {
1535           wYclosest = cl->GetY();
1536           wZclosest = cl->GetZ();
1537
1538           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1539           Double_t h01 = GetTiltFactor(cl);
1540           Int_t det = cl->GetDetector();
1541           Int_t plane = fGeom->GetPlane(det);
1542
1543           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1544           //if(!t.Update(cl,maxChi2,index,h01)) {
1545             if(!tryAgain--) return 0;
1546           }
1547           else tryAgain=fMaxGap;
1548         }
1549         else {
1550           if (tryAgain==0) break;
1551           tryAgain--;
1552
1553           //if (minDY < 1000000 && isNewLayer)
1554             //cout << "\t" << nRefPlane << "\t" << "\t" << t.GetNRotate() <<  "\t" <<
1555             //  road << "\t" << minDY << "\t" << chi2 << "\t" << wChi2 << "\t" << maxChi2 << endl;
1556
1557         }
1558
1559         isNewLayer = kFALSE;
1560
1561         /*
1562         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1563
1564           printf(" %f", wIndex);       //1
1565           printf(" %f", wTB);          //2
1566           printf(" %f", wYrt);         //3
1567           printf(" %f", wYclosest);    //4
1568           printf(" %f", wYcorrect);    //5
1569           printf(" %f", wYwindow);     //6
1570           printf(" %f", wZrt);         //7
1571           printf(" %f", wZclosest);    //8
1572           printf(" %f", wZcorrect);    //9
1573           printf(" %f", wZwindow);     //10
1574           printf(" %f", wPx);          //11
1575           printf(" %f", wPy);          //12
1576           printf(" %f", wPz);          //13
1577           printf(" %f", wSigmaC2*1000000);  //14
1578           printf(" %f", wSigmaTgl2*1000);   //15
1579           printf(" %f", wSigmaY2);     //16
1580           //      printf(" %f", wSigmaZ2);     //17
1581           printf(" %f", wChi2);     //17
1582           printf(" %f", wC);           //18
1583           printf("\n");
1584         }
1585         */
1586       }
1587     }
1588   }
1589   if (nr<outerTB)
1590     t.SetStop(kTRUE);
1591   else
1592     t.SetStop(kFALSE);
1593   return expectedNumberOfClusters;
1594
1595
1596 }
1597
1598 //---------------------------------------------------------------------------
1599 Int_t AliTRDtracker::Refit(AliTRDtrack& t, Int_t rf)
1600 {
1601   // Starting from current position on track=t this function tries
1602   // to extrapolate the track up to timeBin=0 and to reuse already
1603   // assigned clusters. Returns the number of clusters
1604   // expected to be found in sensitive layers
1605   // get indices of assigned clusters for each layer
1606   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
1607
1608   Int_t iCluster[90];
1609   for (Int_t i = 0; i < 90; i++) iCluster[i] = 0;
1610   for (Int_t i = 0; i < t.GetNumberOfClusters(); i++) {
1611     Int_t index = t.GetClusterIndex(i);
1612     AliTRDcluster *cl=(AliTRDcluster*) GetCluster(index);
1613     if (!cl) continue;
1614     Int_t detector=cl->GetDetector();
1615     Int_t localTimeBin=cl->GetLocalTimeBin();
1616     Int_t sector=fGeom->GetSector(detector);
1617     Int_t plane=fGeom->GetPlane(detector);
1618
1619     Int_t trackingSector = CookSectorIndex(sector);
1620
1621     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1622     if(gtb < 0) continue;
1623     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1624     iCluster[layer] = index;
1625   }
1626   t.ResetClusters();
1627
1628   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1629
1630   Double_t alpha=t.GetAlpha();
1631   alpha = TVector2::Phi_0_2pi(alpha);
1632
1633   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1634   Double_t radLength, rho, x, dx, y, ymax, z;
1635
1636   Int_t expectedNumberOfClusters = 0;
1637   Bool_t lookForCluster;
1638
1639   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1640
1641
1642   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
1643
1644     y = t.GetY(); z = t.GetZ();
1645
1646     // first propagate to the inner surface of the current time bin
1647     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1648     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1649     if(!t.PropagateTo(x,radLength,rho)) break;
1650     y = t.GetY();
1651     ymax = x*TMath::Tan(0.5*alpha);
1652     if (y > ymax) {
1653       s = (s+1) % ns;
1654       if (!t.Rotate(alpha)) break;
1655       if(!t.PropagateTo(x,radLength,rho)) break;
1656     } else if (y <-ymax) {
1657       s = (s-1+ns) % ns;
1658       if (!t.Rotate(-alpha)) break;
1659       if(!t.PropagateTo(x,radLength,rho)) break;
1660     }
1661
1662     y = t.GetY(); z = t.GetZ();
1663
1664     // now propagate to the middle plane of the next time bin
1665     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1666     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1667     if(!t.PropagateTo(x,radLength,rho)) break;
1668     y = t.GetY();
1669     ymax = x*TMath::Tan(0.5*alpha);
1670     if (y > ymax) {
1671       s = (s+1) % ns;
1672       if (!t.Rotate(alpha)) break;
1673       if(!t.PropagateTo(x,radLength,rho)) break;
1674     } else if (y <-ymax) {
1675       s = (s-1+ns) % ns;
1676       if (!t.Rotate(-alpha)) break;
1677       if(!t.PropagateTo(x,radLength,rho)) break;
1678     }
1679
1680     if(lookForCluster) expectedNumberOfClusters++;
1681
1682     // use assigned cluster
1683     if (!iCluster[nr-1]) continue;
1684     AliTRDcluster *cl=(AliTRDcluster*)GetCluster(iCluster[nr-1]);
1685     Double_t h01 = GetTiltFactor(cl);
1686     Double_t chi2=t.GetPredictedChi2(cl, h01);
1687     t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1688     t.Update(cl,chi2,iCluster[nr-1],h01);
1689   }
1690
1691   return expectedNumberOfClusters;
1692 }
1693
1694 //___________________________________________________________________
1695
1696 Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1697 {
1698   // Starting from current radial position of track <t> this function
1699   // extrapolates the track up to radial position <xToGo>.
1700   // Returns 1 if track reaches the plane, and 0 otherwise
1701
1702   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1703
1704   Double_t alpha=t.GetAlpha();
1705
1706   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1707   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1708
1709   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1710
1711   Bool_t lookForCluster;
1712   Double_t radLength, rho, x, dx, y, ymax, z;
1713
1714   x = t.GetX();
1715
1716   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1717
1718   Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1719
1720   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) {
1721
1722     y = t.GetY(); z = t.GetZ();
1723
1724     // first propagate to the outer surface of the current time bin
1725     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1726     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
1727     if(!t.PropagateTo(x,radLength,rho)) return 0;
1728     y = t.GetY();
1729     ymax = x*TMath::Tan(0.5*alpha);
1730     if (y > ymax) {
1731       s = (s+1) % ns;
1732       if (!t.Rotate(alpha)) return 0;
1733     } else if (y <-ymax) {
1734       s = (s-1+ns) % ns;
1735       if (!t.Rotate(-alpha)) return 0;
1736     }
1737     if(!t.PropagateTo(x,radLength,rho)) return 0;
1738
1739     y = t.GetY(); z = t.GetZ();
1740
1741     // now propagate to the middle plane of the next time bin
1742     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1743     x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
1744     if(!t.PropagateTo(x,radLength,rho)) return 0;
1745     y = t.GetY();
1746     ymax = x*TMath::Tan(0.5*alpha);
1747     if (y > ymax) {
1748       s = (s+1) % ns;
1749       if (!t.Rotate(alpha)) return 0;
1750     } else if (y <-ymax) {
1751       s = (s-1+ns) % ns;
1752       if (!t.Rotate(-alpha)) return 0;
1753     }
1754     if(!t.PropagateTo(x,radLength,rho)) return 0;
1755   }
1756   return 1;
1757 }
1758
1759 //___________________________________________________________________
1760
1761 Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1762 {
1763   // Starting from current radial position of track <t> this function
1764   // extrapolates the track up to radial position of the outermost
1765   // padrow of the TPC.
1766   // Returns 1 if track reaches the TPC, and 0 otherwise
1767
1768   //Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1769
1770   Double_t alpha=t.GetAlpha();
1771   alpha = TVector2::Phi_0_2pi(alpha);
1772
1773   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1774
1775   Bool_t lookForCluster;
1776   Double_t radLength, rho, x, dx, y, /*ymax,*/ z;
1777
1778   x = t.GetX();
1779
1780   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1781   Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1782
1783   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) {
1784
1785     y = t.GetY();
1786     z = t.GetZ();
1787
1788     // first propagate to the outer surface of the current time bin
1789     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1790     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2;
1791
1792     if(!t.PropagateTo(x,radLength,rho)) return 0;
1793     AdjustSector(&t);
1794     if(!t.PropagateTo(x,radLength,rho)) return 0;
1795
1796     y = t.GetY();
1797     z = t.GetZ();
1798
1799     // now propagate to the middle plane of the next time bin
1800     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1801     x = fTrSec[s]->GetLayer(nr-1)->GetX();
1802
1803     if(!t.PropagateTo(x,radLength,rho)) return 0;
1804     AdjustSector(&t);
1805     if(!t.PropagateTo(x,radLength,rho)) return 0;
1806   }
1807   return 1;
1808 }
1809
1810 void AliTRDtracker::LoadEvent()
1811 {
1812   // Fills clusters into TRD tracking_sectors
1813   // Note that the numbering scheme for the TRD tracking_sectors
1814   // differs from that of TRD sectors
1815
1816   ReadClusters(fClusters);
1817   Int_t ncl=fClusters->GetEntriesFast();
1818   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1819
1820   UInt_t index;
1821   while (ncl--) {
1822 //    printf("\r %d left  ",ncl);
1823     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1824     Int_t detector=c->GetDetector();
1825     Int_t localTimeBin=c->GetLocalTimeBin();
1826     Int_t sector=fGeom->GetSector(detector);
1827     Int_t plane=fGeom->GetPlane(detector);
1828
1829     Int_t trackingSector = CookSectorIndex(sector);
1830
1831     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1832     if(gtb < 0) continue;
1833     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1834
1835     index=ncl;
1836     fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
1837   }
1838   printf("\r\n");
1839
1840 }
1841
1842 //_____________________________________________________________________________
1843 Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1844 {
1845   // Fills clusters into TRD tracking_sectors
1846   // Note that the numbering scheme for the TRD tracking_sectors
1847   // differs from that of TRD sectors
1848
1849   if (ReadClusters(fClusters,cTree)) {
1850      Error("LoadClusters","Problem with reading the clusters !");
1851      return 1;
1852   }
1853   Int_t ncl=fClusters->GetEntriesFast();
1854   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1855
1856   UInt_t index;
1857   for (Int_t ichamber=0;ichamber<5;ichamber++)
1858     for (Int_t isector=0;isector<18;isector++){
1859       fHoles[ichamber][isector]=kTRUE;
1860     }
1861
1862
1863   while (ncl--) {
1864 //    printf("\r %d left  ",ncl);
1865     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1866     Int_t detector=c->GetDetector();
1867     Int_t localTimeBin=c->GetLocalTimeBin();
1868     Int_t sector=fGeom->GetSector(detector);
1869     Int_t plane=fGeom->GetPlane(detector);
1870
1871     Int_t trackingSector = CookSectorIndex(sector);
1872     if (c->GetLabel(0)>0){
1873       Int_t chamber = fGeom->GetChamber(detector);
1874       fHoles[chamber][trackingSector]=kFALSE;
1875     }
1876
1877     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1878     if(gtb < 0) continue;
1879     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1880
1881     index=ncl;
1882     fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
1883   }
1884   printf("\r\n");
1885   //
1886   //
1887   /*
1888   for (Int_t isector=0;isector<18;isector++){
1889     for (Int_t ichamber=0;ichamber<5;ichamber++)
1890       if (fHoles[ichamber][isector]!=fGeom->IsHole(0,ichamber,17-isector))
1891         printf("Problem \t%d\t%d\t%d\t%d\n",isector,ichamber,fHoles[ichamber][isector],
1892              fGeom->IsHole(0,ichamber,17-isector));
1893   }
1894   */
1895   return 0;
1896 }
1897
1898 //_____________________________________________________________________________
1899 void AliTRDtracker::UnloadEvent()
1900 {
1901   //
1902   // Clears the arrays of clusters and tracks. Resets sectors and timebins
1903   //
1904
1905   Int_t i, nentr;
1906
1907   nentr = fClusters->GetEntriesFast();
1908   for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1909
1910   nentr = fSeeds->GetEntriesFast();
1911   for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1912
1913   nentr = fTracks->GetEntriesFast();
1914   for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1915
1916   Int_t nsec = AliTRDgeometry::kNsect;
1917
1918   for (i = 0; i < nsec; i++) {
1919     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1920       fTrSec[i]->GetLayer(pl)->Clear();
1921     }
1922   }
1923
1924 }
1925
1926 //__________________________________________________________________________
1927 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1928 {
1929   // Creates track seeds using clusters in timeBins=i1,i2
1930
1931   if(turn > 2) {
1932     cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
1933     return;
1934   }
1935
1936   Double_t x[5], c[15];
1937   Int_t maxSec=AliTRDgeometry::kNsect;
1938
1939   Double_t alpha=AliTRDgeometry::GetAlpha();
1940   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1941   Double_t cs=cos(alpha), sn=sin(alpha);
1942   Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1943
1944
1945   Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1946   Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1947
1948   Double_t x1 =fTrSec[0]->GetX(i1);
1949   Double_t xx2=fTrSec[0]->GetX(i2);
1950
1951   for (Int_t ns=0; ns<maxSec; ns++) {
1952
1953     Int_t nl2 = *(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
1954     Int_t nl=(*fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
1955     Int_t nm=(*fTrSec[ns]->GetLayer(i2));
1956     Int_t nu=(*fTrSec[(ns+1)%maxSec]->GetLayer(i2));
1957     Int_t nu2=(*fTrSec[(ns+2)%maxSec]->GetLayer(i2));
1958
1959     AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1960
1961     for (Int_t is=0; is < r1; is++) {
1962       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
1963
1964       for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
1965
1966         const AliTRDcluster *cl;
1967         Double_t x2,   y2,   z2;
1968         Double_t x3=0., y3=0.;
1969
1970         if (js<nl2) {
1971           if(turn != 2) continue;
1972           AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
1973           cl=r2[js];
1974           y2=cl->GetY(); z2=cl->GetZ();
1975
1976           x2= xx2*cs2+y2*sn2;
1977           y2=-xx2*sn2+y2*cs2;
1978         }
1979         else if (js<nl2+nl) {
1980           if(turn != 1) continue;
1981           AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
1982           cl=r2[js-nl2];
1983           y2=cl->GetY(); z2=cl->GetZ();
1984
1985           x2= xx2*cs+y2*sn;
1986           y2=-xx2*sn+y2*cs;
1987         }
1988         else if (js<nl2+nl+nm) {
1989           if(turn != 1) continue;
1990           AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1991           cl=r2[js-nl2-nl];
1992           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1993         }
1994         else if (js<nl2+nl+nm+nu) {
1995           if(turn != 1) continue;
1996           AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%maxSec]->GetLayer(i2));
1997           cl=r2[js-nl2-nl-nm];
1998           y2=cl->GetY(); z2=cl->GetZ();
1999
2000           x2=xx2*cs-y2*sn;
2001           y2=xx2*sn+y2*cs;
2002         }
2003         else {
2004           if(turn != 2) continue;
2005           AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%maxSec]->GetLayer(i2));
2006           cl=r2[js-nl2-nl-nm-nu];
2007           y2=cl->GetY(); z2=cl->GetZ();
2008
2009           x2=xx2*cs2-y2*sn2;
2010           y2=xx2*sn2+y2*cs2;
2011         }
2012
2013         if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue;
2014
2015         Double_t zz=z1 - z1/x1*(x1-x2);
2016
2017         if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue;
2018
2019         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
2020         if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
2021
2022         x[0]=y1;
2023         x[1]=z1;
2024         x[4]=f1trd(x1,y1,x2,y2,x3,y3);
2025
2026         if (TMath::Abs(x[4]) > fgkMaxSeedC) continue;
2027
2028         x[2]=f2trd(x1,y1,x2,y2,x3,y3);
2029
2030         if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
2031
2032         x[3]=f3trd(x1,y1,x2,y2,z1,z2);
2033
2034         if (TMath::Abs(x[3]) > fgkMaxSeedTan) continue;
2035
2036         Double_t a=asin(x[2]);
2037         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
2038
2039         if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue;
2040
2041         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
2042         Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2();
2043         Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ;
2044
2045         // Tilt changes
2046         Double_t h01 = GetTiltFactor(r1[is]);
2047         Double_t xuFactor = 100.;
2048         if(fNoTilt) {
2049           h01 = 0;
2050           xuFactor = 1;
2051         }
2052
2053         sy1=sy1+sz1*h01*h01;
2054         Double_t syz=sz1*(-h01);
2055         // end of tilt changes
2056
2057         Double_t f40=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2058         Double_t f42=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2059         Double_t f43=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2060         Double_t f20=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2061         Double_t f22=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2062         Double_t f23=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2063         Double_t f30=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2064         Double_t f31=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2065         Double_t f32=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2066         Double_t f34=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2067
2068
2069         c[0]=sy1;
2070         //        c[1]=0.;       c[2]=sz1;
2071         c[1]=syz;       c[2]=sz1*xuFactor;
2072         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2073         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
2074                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2075         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2076         c[13]=f30*sy1*f40+f32*sy2*f42;
2077         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2078
2079         UInt_t index=r1.GetIndex(is);
2080
2081         AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
2082
2083         Int_t rc=FollowProlongation(*track, i2);
2084
2085         if ((rc < 1) ||
2086             (track->GetNumberOfClusters() <
2087              (outer-inner)*fgkMinClustersInSeed)) delete track;
2088         else {
2089           fSeeds->AddLast(track); fNseeds++;
2090 //          cerr<<"\r found seed "<<fNseeds;
2091         }
2092       }
2093     }
2094   }
2095 }
2096
2097 //_____________________________________________________________________________
2098 Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree)
2099 {
2100   //
2101   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2102   // from the file. The names of the cluster tree and branches
2103   // should match the ones used in AliTRDclusterizer::WriteClusters()
2104   //
2105   TObjArray *clusterArray = new TObjArray(400);
2106
2107   TBranch *branch=ClusterTree->GetBranch("TRDcluster");
2108   if (!branch) {
2109     Error("ReadClusters","Can't get the branch !");
2110     return 1;
2111   }
2112   branch->SetAddress(&clusterArray);
2113
2114   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
2115   printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
2116
2117   // Loop through all entries in the tree
2118   Int_t nbytes;
2119   AliTRDcluster *c = 0;
2120   printf("\n");
2121
2122   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2123
2124     // Import the tree
2125     nbytes += ClusterTree->GetEvent(iEntry);
2126
2127     // Get the number of points in the detector
2128     Int_t nCluster = clusterArray->GetEntriesFast();
2129 //    printf("\r Read %d clusters from entry %d", nCluster, iEntry);
2130
2131     // Loop through all TRD digits
2132     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2133       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2134       AliTRDcluster *co = new AliTRDcluster(*c);
2135       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
2136       Int_t ltb = co->GetLocalTimeBin();
2137       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
2138       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
2139       array->AddLast(co);
2140       delete clusterArray->RemoveAt(iCluster);
2141     }
2142   }
2143
2144   delete clusterArray;
2145
2146   return 0;
2147 }
2148
2149 //______________________________________________________________________
2150 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
2151 {
2152   //
2153   // Reads AliTRDclusters from file <filename>. The names of the cluster
2154   // tree and branches should match the ones used in
2155   // AliTRDclusterizer::WriteClusters()
2156   // if <array> == 0, clusters are added into AliTRDtracker fCluster array
2157   //
2158
2159   TDirectory *savedir=gDirectory;
2160
2161   TFile *file = TFile::Open(filename);
2162   if (!file->IsOpen()) {
2163     cerr<<"Can't open file with TRD clusters"<<endl;
2164     return;
2165   }
2166
2167   Char_t treeName[12];
2168   sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
2169   TTree *clusterTree = (TTree*) gDirectory->Get(treeName);
2170
2171   if (!clusterTree) {
2172      cerr<<"AliTRDtracker::ReadClusters(): ";
2173      cerr<<"can't get a tree with clusters !\n";
2174      return;
2175   }
2176
2177   TObjArray *clusterArray = new TObjArray(400);
2178
2179   clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
2180
2181   Int_t nEntries = (Int_t) clusterTree->GetEntries();
2182   cout<<"found "<<nEntries<<" in clusterTree"<<endl;
2183
2184   // Loop through all entries in the tree
2185   Int_t nbytes;
2186   AliTRDcluster *c = 0;
2187
2188   printf("\n");
2189
2190   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2191
2192     // Import the tree
2193     nbytes += clusterTree->GetEvent(iEntry);
2194
2195     // Get the number of points in the detector
2196     Int_t nCluster = clusterArray->GetEntriesFast();
2197     printf("\n Read %d clusters from entry %d", nCluster, iEntry);
2198
2199     // Loop through all TRD digits
2200     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2201       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2202       AliTRDcluster *co = new AliTRDcluster(*c);
2203       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
2204       Int_t ltb = co->GetLocalTimeBin();
2205       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
2206       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
2207       array->AddLast(co);
2208       delete clusterArray->RemoveAt(iCluster);
2209     }
2210   }
2211
2212   file->Close();
2213   delete clusterArray;
2214   savedir->cd();
2215
2216 }
2217
2218 void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp)
2219 {
2220   //
2221   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2222   // from the file. The names of the cluster tree and branches
2223   // should match the ones used in AliTRDclusterizer::WriteClusters()
2224   //
2225
2226   TDirectory *savedir=gDirectory;
2227
2228   if (inp) {
2229      TFile *in=(TFile*)inp;
2230      if (!in->IsOpen()) {
2231         cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n";
2232         return;
2233      }
2234      else{
2235        in->cd();
2236      }
2237   }
2238
2239   Char_t treeName[12];
2240   sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
2241   TTree *clusterTree = (TTree*) gDirectory->Get(treeName);
2242
2243   TObjArray *clusterArray = new TObjArray(400);
2244
2245   clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
2246
2247   Int_t nEntries = (Int_t) clusterTree->GetEntries();
2248   printf("found %d entries in %s.\n",nEntries,clusterTree->GetName());
2249
2250   // Loop through all entries in the tree
2251   Int_t nbytes;
2252   AliTRDcluster *c = 0;
2253   printf("\n");
2254
2255   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2256
2257     // Import the tree
2258     nbytes += clusterTree->GetEvent(iEntry);
2259
2260     // Get the number of points in the detector
2261     Int_t nCluster = clusterArray->GetEntriesFast();
2262 //    printf("\r Read %d clusters from entry %d", nCluster, iEntry);
2263
2264     // Loop through all TRD digits
2265     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2266       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2267       AliTRDcluster *co = new AliTRDcluster(*c);
2268       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
2269       Int_t ltb = co->GetLocalTimeBin();
2270       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
2271       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
2272       array->AddLast(co);
2273       delete clusterArray->RemoveAt(iCluster);
2274     }
2275   }
2276
2277   delete clusterArray;
2278   savedir->cd();
2279
2280 }
2281
2282 //__________________________________________________________________
2283 void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const
2284 {
2285   //
2286   // This cooks a label. Mmmmh, smells good...
2287   //
2288
2289   Int_t label=123456789, index, i, j;
2290   Int_t ncl=pt->GetNumberOfClusters();
2291   const Int_t kRange = fTrSec[0]->GetOuterTimeBin()+1;
2292
2293   Bool_t labelAdded;
2294
2295   //  Int_t s[kRange][2];
2296   Int_t **s = new Int_t* [kRange];
2297   for (i=0; i<kRange; i++) {
2298     s[i] = new Int_t[2];
2299   }
2300   for (i=0; i<kRange; i++) {
2301     s[i][0]=-1;
2302     s[i][1]=0;
2303   }
2304
2305   Int_t t0,t1,t2;
2306   for (i=0; i<ncl; i++) {
2307     index=pt->GetClusterIndex(i);
2308     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2309     t0=c->GetLabel(0);
2310     t1=c->GetLabel(1);
2311     t2=c->GetLabel(2);
2312   }
2313
2314   for (i=0; i<ncl; i++) {
2315     index=pt->GetClusterIndex(i);
2316     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2317     for (Int_t k=0; k<3; k++) {
2318       label=c->GetLabel(k);
2319       labelAdded=kFALSE; j=0;
2320       if (label >= 0) {
2321         while ( (!labelAdded) && ( j < kRange ) ) {
2322           if (s[j][0]==label || s[j][1]==0) {
2323             s[j][0]=label;
2324             s[j][1]=s[j][1]+1;
2325             labelAdded=kTRUE;
2326           }
2327           j++;
2328         }
2329       }
2330     }
2331   }
2332
2333   Int_t max=0;
2334   label = -123456789;
2335
2336   for (i=0; i<kRange; i++) {
2337     if (s[i][1]>max) {
2338       max=s[i][1]; label=s[i][0];
2339     }
2340   }
2341
2342   for (i=0; i<kRange; i++) {
2343     delete []s[i];
2344   }
2345
2346   delete []s;
2347
2348   if ((1.- Float_t(max)/ncl) > wrong) label=-label;
2349
2350   pt->SetLabel(label);
2351
2352 }
2353
2354
2355 //__________________________________________________________________
2356 void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const
2357 {
2358   //
2359   // Use clusters, but don't abuse them!
2360   //
2361
2362   Int_t ncl=t->GetNumberOfClusters();
2363   for (Int_t i=from; i<ncl; i++) {
2364     Int_t index = t->GetClusterIndex(i);
2365     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2366     c->Use();
2367   }
2368 }
2369
2370
2371 //_____________________________________________________________________
2372 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
2373 {
2374   // Parametrised "expected" error of the cluster reconstruction in Y
2375
2376   Double_t s = 0.08 * 0.08;
2377   return s;
2378 }
2379
2380 //_____________________________________________________________________
2381 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
2382 {
2383   // Parametrised "expected" error of the cluster reconstruction in Z
2384
2385   Double_t s = 9 * 9 /12.;
2386   return s;
2387 }
2388
2389 //_____________________________________________________________________
2390 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
2391 {
2392   //
2393   // Returns radial position which corresponds to time bin <localTB>
2394   // in tracking sector <sector> and plane <plane>
2395   //
2396
2397   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
2398   Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2399   return fTrSec[sector]->GetLayer(pl)->GetX();
2400
2401 }
2402
2403
2404 //_______________________________________________________
2405 AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
2406                Double_t dx, Double_t rho, Double_t radLength, Int_t tbIndex)
2407 {
2408   //
2409   // AliTRDpropagationLayer constructor
2410   //
2411
2412   fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = radLength;
2413   fClusters = NULL; fIndex = NULL; fTimeBinIndex = tbIndex;
2414
2415
2416   for(Int_t i=0; i < (Int_t) kZones; i++) {
2417     fZc[i]=0; fZmax[i] = 0;
2418   }
2419
2420   fYmax = 0;
2421
2422   if(fTimeBinIndex >= 0) {
2423     fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
2424     fIndex = new UInt_t[kMaxClusterPerTimeBin];
2425   }
2426
2427   for (Int_t i=0;i<5;i++) fIsHole[i] = kFALSE;
2428   fHole = kFALSE;
2429   fHoleZc = 0;
2430   fHoleZmax = 0;
2431   fHoleYc = 0;
2432   fHoleYmax = 0;
2433   fHoleRho = 0;
2434   fHoleX0 = 0;
2435
2436 }
2437
2438 //_______________________________________________________
2439 void AliTRDtracker::AliTRDpropagationLayer::SetHole(
2440           Double_t Zmax, Double_t Ymax, Double_t rho,
2441           Double_t radLength, Double_t Yc, Double_t Zc)
2442 {
2443   //
2444   // Sets hole in the layer
2445   //
2446   fHole = kTRUE;
2447   fHoleZc = Zc;
2448   fHoleZmax = Zmax;
2449   fHoleYc = Yc;
2450   fHoleYmax = Ymax;
2451   fHoleRho = rho;
2452   fHoleX0 = radLength;
2453 }
2454
2455
2456 //_______________________________________________________
2457 AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
2458 {
2459   //
2460   // AliTRDtrackingSector Constructor
2461   //
2462
2463   fGeom = geo;
2464   fPar = par;
2465   fGeomSector = gs;
2466   fTzeroShift = 0.13;
2467   fN = 0;
2468   //
2469   // get holes description from geometry
2470   Bool_t holes[AliTRDgeometry::kNcham];
2471   //printf("sector\t%d\t",gs);
2472   for (Int_t icham=0; icham<AliTRDgeometry::kNcham;icham++){
2473     holes[icham] = fGeom->IsHole(0,icham,gs);
2474     //printf("%d",holes[icham]);
2475   }
2476   //printf("\n");
2477
2478   for(UInt_t i=0; i < kMaxTimeBinIndex; i++) fTimeBinIndex[i] = -1;
2479
2480
2481   AliTRDpropagationLayer* ppl;
2482
2483   Double_t x, xin, xout, dx, rho, radLength;
2484   Int_t    steps;
2485
2486   // set time bins in the gas of the TPC
2487
2488   xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
2489   rho = 0.9e-3;  radLength = 28.94;
2490
2491   for(Int_t i=0; i<steps; i++) {
2492     x = xin + i*dx + dx/2;
2493     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2494     InsertLayer(ppl);
2495   }
2496
2497   // set time bins in the outer field cage vessel
2498
2499   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2500   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2501   InsertLayer(ppl);
2502
2503   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2504   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2505   InsertLayer(ppl);
2506
2507   dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
2508   steps = 5; dx = (xout - xin)/steps;
2509   for(Int_t i=0; i<steps; i++) {
2510     x = xin + i*dx + dx/2;
2511     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2512     InsertLayer(ppl);
2513   }
2514
2515   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2516   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2517   InsertLayer(ppl);
2518
2519   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2520   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2521   InsertLayer(ppl);
2522
2523
2524   // set time bins in CO2
2525
2526   xin = xout; xout = 275.0;
2527   steps = 50; dx = (xout - xin)/steps;
2528   rho = 1.977e-3;  radLength = 36.2;
2529
2530   for(Int_t i=0; i<steps; i++) {
2531     x = xin + i*dx + dx/2;
2532     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2533     InsertLayer(ppl);
2534   }
2535
2536   // set time bins in the outer containment vessel
2537
2538   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2539   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2540   InsertLayer(ppl);
2541
2542   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2543   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2544   InsertLayer(ppl);
2545
2546   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2547   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2548   InsertLayer(ppl);
2549
2550   dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
2551   steps = 10; dx = (xout - xin)/steps;
2552   for(Int_t i=0; i<steps; i++) {
2553     x = xin + i*dx + dx/2;
2554     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2555     InsertLayer(ppl);
2556   }
2557
2558   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2559   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2560   InsertLayer(ppl);
2561
2562   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2563   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2564   InsertLayer(ppl);
2565
2566   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2567   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2568   InsertLayer(ppl);
2569
2570   Double_t xtrd = (Double_t) fGeom->Rmin();
2571
2572   // add layers between TPC and TRD (Air temporarily)
2573   xin = xout; xout = xtrd;
2574   steps = 50; dx = (xout - xin)/steps;
2575   rho = 1.2e-3;  radLength = 36.66;
2576
2577   for(Int_t i=0; i<steps; i++) {
2578     x = xin + i*dx + dx/2;
2579     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2580     InsertLayer(ppl);
2581   }
2582
2583
2584   //  Double_t alpha=AliTRDgeometry::GetAlpha();
2585
2586   // add layers for each of the planes
2587
2588   Double_t dxRo = (Double_t) fGeom->CroHght();    // Rohacell
2589   Double_t dxSpace = (Double_t) fGeom->Cspace();  // Spacing between planes
2590   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2591   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2592   Double_t dxRad = (Double_t) fGeom->CraHght();   // Radiator
2593   Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo;
2594   Double_t dxPlane = dxTEC + dxSpace;
2595
2596   Int_t tb, tbIndex;
2597   const Int_t  kNchambers = AliTRDgeometry::Ncham();
2598   Double_t  ymax = 0;
2599   //, holeYmax = 0;
2600   Double_t ymaxsensitive=0;
2601   Double_t *zc = new Double_t[kNchambers];
2602   Double_t *zmax = new Double_t[kNchambers];
2603   Double_t *zmaxsensitive = new Double_t[kNchambers];
2604   //  Double_t  holeZmax = 1000.;   // the whole sector is missing
2605
2606   for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2607     //
2608     // Radiator
2609     xin = xtrd + plane * dxPlane; xout = xin + dxRad;
2610     steps = 12; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
2611     for(Int_t i=0; i<steps; i++) {
2612       x = xin + i*dx + dx/2;
2613       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2614       InsertLayer(ppl);
2615     }
2616
2617     ymax          = fGeom->GetChamberWidth(plane)/2.;
2618     ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
2619
2620     for(Int_t ch = 0; ch < kNchambers; ch++) {
2621       zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
2622       Float_t pad = fPar->GetRowPadSize(plane,ch,0);
2623       Float_t row0 = fPar->GetRow0(plane,ch,0);
2624       Int_t nPads = fPar->GetRowMax(plane,ch,0);
2625       zmaxsensitive[ch] = Float_t(nPads)*pad/2.;
2626       //      zc[ch] = (pad * nPads)/2 + row0 - pad/2;
2627       zc[ch] = (pad * nPads)/2 + row0;
2628       //zc[ch] = row0+zmax[ch]-AliTRDgeometry::RpadW();
2629
2630     }
2631
2632     dx = fPar->GetTimeBinSize();
2633     rho = 0.00295 * 0.85; radLength = 11.0;
2634
2635     Double_t x0 = (Double_t) fPar->GetTime0(plane);
2636     Double_t xbottom = x0 - dxDrift;
2637     Double_t xtop = x0 + dxAmp;
2638     //
2639     // Amplification region
2640     steps = (Int_t) (dxAmp/dx);
2641
2642     for(tb = 0; tb < steps; tb++) {
2643       x = x0 + tb * dx + dx/2;
2644       tbIndex = CookTimeBinIndex(plane, -tb-1);
2645       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2646       ppl->SetYmax(ymax,ymaxsensitive);
2647       ppl->SetZ(zc, zmax, zmaxsensitive);
2648       ppl->SetHoles(holes);
2649       InsertLayer(ppl);
2650     }
2651     tbIndex = CookTimeBinIndex(plane, -steps);
2652     x = (x + dx/2 + xtop)/2;
2653     dx = 2*(xtop-x);
2654     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2655     ppl->SetYmax(ymax,ymaxsensitive);
2656     ppl->SetZ(zc, zmax,zmaxsensitive);
2657     ppl->SetHoles(holes);
2658     InsertLayer(ppl);
2659
2660     // Drift region
2661     dx = fPar->GetTimeBinSize();
2662     steps = (Int_t) (dxDrift/dx);
2663
2664     for(tb = 0; tb < steps; tb++) {
2665       x = x0 - tb * dx - dx/2;
2666       tbIndex = CookTimeBinIndex(plane, tb);
2667
2668       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2669       ppl->SetYmax(ymax,ymaxsensitive);
2670       ppl->SetZ(zc, zmax, zmaxsensitive);
2671       ppl->SetHoles(holes);
2672       InsertLayer(ppl);
2673     }
2674     tbIndex = CookTimeBinIndex(plane, steps);
2675     x = (x - dx/2 + xbottom)/2;
2676     dx = 2*(x-xbottom);
2677     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2678     ppl->SetYmax(ymax,ymaxsensitive);
2679     ppl->SetZ(zc, zmax, zmaxsensitive);
2680     ppl->SetHoles(holes);
2681     InsertLayer(ppl);
2682
2683     // Pad Plane
2684     xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; radLength = 33.0;
2685     ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2686     ppl->SetYmax(ymax,ymaxsensitive);
2687     ppl->SetZ(zc, zmax,zmax);
2688     ppl->SetHoles(holes);
2689     InsertLayer(ppl);
2690
2691     // Rohacell
2692     xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
2693     steps = 5; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
2694     for(Int_t i=0; i<steps; i++) {
2695       x = xin + i*dx + dx/2;
2696       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2697       ppl->SetYmax(ymax,ymaxsensitive);
2698       ppl->SetZ(zc, zmax,zmax);
2699       ppl->SetHoles(holes);
2700       InsertLayer(ppl);
2701     }
2702
2703     // Space between the chambers, air
2704     xin = xout; xout = xtrd + (plane + 1) * dxPlane;
2705     steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
2706     for(Int_t i=0; i<steps; i++) {
2707       x = xin + i*dx + dx/2;
2708       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2709       InsertLayer(ppl);
2710     }
2711   }
2712
2713   // Space between the TRD and RICH
2714   Double_t xRICH = 500.;
2715   xin = xout; xout = xRICH;
2716   steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
2717   for(Int_t i=0; i<steps; i++) {
2718     x = xin + i*dx + dx/2;
2719     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2720     InsertLayer(ppl);
2721   }
2722
2723   MapTimeBinLayers();
2724   delete [] zc;
2725   delete [] zmax;
2726
2727 }
2728
2729 //______________________________________________________
2730
2731 Int_t  AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t localTB) const
2732 {
2733   //
2734   // depending on the digitization parameters calculates "global"
2735   // time bin index for timebin <localTB> in plane <plane>
2736   //
2737
2738   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2739   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2740   Double_t dx = (Double_t) fPar->GetTimeBinSize();
2741
2742   Int_t tbAmp = fPar->GetTimeBefore();
2743   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
2744   if(kTRUE) maxAmp = 0;   // intentional until we change parameter class
2745   Int_t tbDrift = fPar->GetTimeMax();
2746   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
2747
2748   Int_t tbPerPlane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
2749
2750   Int_t gtb = (plane+1) * tbPerPlane - localTB - 1 - TMath::Min(tbAmp,maxAmp);
2751
2752   if((localTB < 0) &&
2753      (TMath::Abs(localTB) > TMath::Min(tbAmp,maxAmp))) return -1;
2754   if(localTB >= TMath::Min(tbDrift,maxDrift)) return -1;
2755
2756   return gtb;
2757
2758
2759 }
2760
2761 //______________________________________________________
2762
2763 void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
2764 {
2765   //
2766   // For all sensitive time bins sets corresponding layer index
2767   // in the array fTimeBins
2768   //
2769
2770   Int_t index;
2771
2772   for(Int_t i = 0; i < fN; i++) {
2773     index = fLayers[i]->GetTimeBinIndex();
2774
2775     //    printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2776
2777     if(index < 0) continue;
2778     if(index >= (Int_t) kMaxTimeBinIndex) {
2779       printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2780       printf("    index %d exceeds allowed maximum of %d!\n",
2781              index, kMaxTimeBinIndex-1);
2782       continue;
2783     }
2784     fTimeBinIndex[index] = i;
2785   }
2786
2787   Double_t x1, dx1, x2, dx2, gap;
2788
2789   for(Int_t i = 0; i < fN-1; i++) {
2790     x1 = fLayers[i]->GetX();
2791     dx1 = fLayers[i]->GetdX();
2792     x2 = fLayers[i+1]->GetX();
2793     dx2 = fLayers[i+1]->GetdX();
2794     gap = (x2 - dx2/2) - (x1 + dx1/2);
2795     if(gap < -0.01) {
2796       printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2797       printf("             %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2798     }
2799     if(gap > 0.01) {
2800       printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2801       printf("             (%f - %f) - (%f + %f) = %f\n",
2802              x2, dx2/2, x1, dx1, gap);
2803     }
2804   }
2805 }
2806
2807
2808 //______________________________________________________
2809
2810
2811 Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2812 {
2813   //
2814   // Returns the number of time bin which in radial position is closest to <x>
2815   //
2816
2817   if(x >= fLayers[fN-1]->GetX()) return fN-1;
2818   if(x <= fLayers[0]->GetX()) return 0;
2819
2820   Int_t b=0, e=fN-1, m=(b+e)/2;
2821   for (; b<e; m=(b+e)/2) {
2822     if (x > fLayers[m]->GetX()) b=m+1;
2823     else e=m;
2824   }
2825   if(TMath::Abs(x - fLayers[m]->GetX()) >
2826      TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2827   else return m;
2828
2829 }
2830
2831 //______________________________________________________
2832
2833 Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
2834 {
2835   //
2836   // Returns number of the innermost SENSITIVE propagation layer
2837   //
2838
2839   return GetLayerNumber(0);
2840 }
2841
2842 //______________________________________________________
2843
2844 Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
2845 {
2846   //
2847   // Returns number of the outermost SENSITIVE time bin
2848   //
2849
2850   return GetLayerNumber(GetNumberOfTimeBins() - 1);
2851 }
2852
2853 //______________________________________________________
2854
2855 Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
2856 {
2857   //
2858   // Returns number of SENSITIVE time bins
2859   //
2860
2861   Int_t tb, layer;
2862   for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
2863     layer = GetLayerNumber(tb);
2864     if(layer>=0) break;
2865   }
2866   return tb+1;
2867 }
2868
2869 //______________________________________________________
2870
2871 void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2872 {
2873   //
2874   // Insert layer <pl> in fLayers array.
2875   // Layers are sorted according to X coordinate.
2876
2877   if ( fN == ((Int_t) kMaxLayersPerSector)) {
2878     printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2879     return;
2880   }
2881   if (fN==0) {fLayers[fN++] = pl; return;}
2882   Int_t i=Find(pl->GetX());
2883
2884   memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2885   fLayers[i]=pl; fN++;
2886
2887 }
2888
2889 //______________________________________________________
2890
2891 Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
2892 {
2893   //
2894   // Returns index of the propagation layer nearest to X
2895   //
2896
2897   if (x <= fLayers[0]->GetX()) return 0;
2898   if (x > fLayers[fN-1]->GetX()) return fN;
2899   Int_t b=0, e=fN-1, m=(b+e)/2;
2900   for (; b<e; m=(b+e)/2) {
2901     if (x > fLayers[m]->GetX()) b=m+1;
2902     else e=m;
2903   }
2904   return m;
2905 }
2906
2907 //______________________________________________________
2908 void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
2909 {
2910   //
2911   // set centers and the width of sectors
2912   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
2913     fZc[icham] = center[icham];
2914     fZmax[icham] = w[icham];
2915     fZmaxSensitive[icham] = wsensitive[icham];
2916     //   printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
2917   }
2918 }
2919 //______________________________________________________
2920
2921 void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
2922 {
2923   //
2924   // set centers and the width of sectors
2925   fHole = kFALSE;
2926   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
2927     fIsHole[icham] = holes[icham];
2928     if (holes[icham]) fHole = kTRUE;
2929   }
2930 }
2931
2932
2933
2934 void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
2935         Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength,
2936         Bool_t &lookForCluster) const
2937 {
2938   //
2939   // Returns radial step <dx>, density <rho>, rad. length <radLength>,
2940   // and sensitivity <lookForCluster> in point <y,z>
2941   //
2942
2943   dx  = fdX;
2944   rho = fRho;
2945   radLength  = fX0;
2946   lookForCluster = kFALSE;
2947   //
2948   // check dead regions in sensitive volume
2949   if(fTimeBinIndex >= 0) {
2950     //
2951     Int_t zone=-1;
2952     for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
2953       if  (TMath::Abs(z - fZc[ch]) < fZmaxSensitive[ch]){
2954         zone = ch;
2955         lookForCluster = !(fIsHole[zone]);
2956         if(TMath::Abs(y) > fYmaxSensitive){
2957           lookForCluster = kFALSE;
2958         }
2959         if (fIsHole[zone]) {
2960           //if hole
2961           rho = 1.29e-3;
2962           radLength = 36.66;
2963         }
2964       }
2965     }
2966     return;
2967   }
2968   //
2969   //
2970   // check hole
2971   if (fHole==kFALSE) return;
2972   //
2973   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
2974     if  (TMath::Abs(z - fZc[ch]) < fZmax[ch]){
2975       if (fIsHole[ch]) {
2976         //if hole
2977         rho = 1.29e-3;
2978         radLength = 36.66;
2979       }
2980     }
2981   }
2982   return;
2983 }
2984
2985 Int_t  AliTRDtracker::AliTRDpropagationLayer::GetZone( Double_t z) const
2986 {
2987   //
2988   //
2989   if (fTimeBinIndex < 0) return -20;  //unknown
2990   Int_t zone=-10;   // dead zone
2991   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
2992     if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
2993       zone = ch;
2994   }
2995   return zone;
2996 }
2997
2998
2999 //______________________________________________________
3000
3001 void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
3002                                                           UInt_t index) {
3003
3004 // Insert cluster in cluster array.
3005 // Clusters are sorted according to Y coordinate.
3006
3007   if(fTimeBinIndex < 0) {
3008     printf("*** attempt to insert cluster into non-sensitive time bin!\n");
3009     return;
3010   }
3011
3012   if (fN== (Int_t) kMaxClusterPerTimeBin) {
3013     printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
3014     return;
3015   }
3016   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
3017   Int_t i=Find(c->GetY());
3018   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3019   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
3020   fIndex[i]=index; fClusters[i]=c; fN++;
3021 }
3022
3023 //______________________________________________________
3024
3025 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
3026
3027 // Returns index of the cluster nearest in Y
3028
3029   if (y <= fClusters[0]->GetY()) return 0;
3030   if (y > fClusters[fN-1]->GetY()) return fN;
3031   Int_t b=0, e=fN-1, m=(b+e)/2;
3032   for (; b<e; m=(b+e)/2) {
3033     if (y > fClusters[m]->GetY()) b=m+1;
3034     else e=m;
3035   }
3036   return m;
3037 }
3038
3039 //---------------------------------------------------------
3040
3041 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
3042 //
3043 //  Returns correction factor for tilted pads geometry
3044 //
3045
3046   Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
3047   Int_t det = c->GetDetector();
3048   Int_t plane = fGeom->GetPlane(det);
3049
3050   if((plane == 1) || (plane == 3) || (plane == 5)) h01=-h01;
3051
3052   if(fNoTilt) h01 = 0;
3053
3054   return h01;
3055 }
3056
3057
3058
3059
3060
3061