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