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