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