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