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 //  Based on Kalman filltering approach                                      //
  22 //                                                                           //
  23 //  Authors:                                                                 //
  24 //     M. Ivanov  (Marian.Ivanov@cern.ch)                                    //
  25 //     Y. Belikov (Jouri.Belikov@cern.ch)                                    //
  26 //                                                                           //
  27 ///////////////////////////////////////////////////////////////////////////////
  28
  29 #include <Riostream.h>
  30 #include <TFile.h>
  31 #include <TBranch.h>
  32 #include <TTree.h>
  33 #include <TObjArray.h>
  34 #include <TTreeStream.h>
  35 #include <TGraph.h>
  36 #include <TLinearFitter.h>
  37
  38 #include "AliESD.h"
  39 #include "AliAlignObj.h"
  40 #include "AliRieman.h"
  41 #include "AliTrackPointArray.h"
  42
  43 #include "AliTRDgeometry.h"
  44 #include "AliTRDpadPlane.h"
  45 #include "AliTRDgeometry.h"
  46 #include "AliTRDcluster.h"
  47 #include "AliTRDtrack.h"
  48 #include "AliTRDseed.h"
  49 #include "AliTRDcalibDB.h"
  50 #include "AliTRDCommonParam.h"
  51 #include "AliTRDtracker.h"
  52 #include "AliTRDReconstructor.h"
  53 #include "AliTRDCalibra.h"
  54
  55 ClassImp(AliTRDtracker)
  56
  57   const  Float_t  AliTRDtracker::fgkMinClustersInTrack =  0.5;
  58   const  Float_t  AliTRDtracker::fgkLabelFraction      =  0.8;
  59   const  Double_t AliTRDtracker::fgkMaxChi2            = 12.0;
  60   const  Double_t AliTRDtracker::fgkMaxSnp             =  0.95; // Corresponds to tan = 3
  61   const  Double_t AliTRDtracker::fgkMaxStep            =  2.0;  // Maximal step size in propagation
  62
  63 //_____________________________________________________________________________
  64 AliTRDtracker::AliTRDtracker()
  65   :AliTracker()
  66   ,fGeom(0)
  67   ,fNclusters(0)
  68   ,fClusters(0)
  69   ,fNseeds(0)
  70   ,fSeeds(0)
  71   ,fNtracks(0)
  72   ,fTracks(0)
  73   ,fTimeBinsPerPlane(0)
  74   ,fAddTRDseeds(kFALSE)
  75   ,fNoTilt(kFALSE)
  76   ,fDebugStreamer(0)
  77 {
  78   //
  79   // Default constructor
  80   //
  81
  82   for (Int_t i = 0; i < kTrackingSectors; i++) {
  83     fTrSec[i] = 0;
  84   }
  85   for (Int_t j = 0; j <  5; j++) {
  86     for (Int_t k = 0; k < 18; k++) {
  87       fHoles[j][k] = kFALSE;
  88     }
  89   }
  90
  91 }
  92
  93 //_____________________________________________________________________________
  94 AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
  95   :AliTracker(t)
  96   ,fGeom(0)
  97   ,fNclusters(0)
  98   ,fClusters(0)
  99   ,fNseeds(0)
 100   ,fSeeds(0)
 101   ,fNtracks(0)
 102   ,fTracks(0)
 103   ,fTimeBinsPerPlane(0)
 104   ,fAddTRDseeds(kFALSE)
 105   ,fNoTilt(kFALSE)
 106   ,fDebugStreamer(0)
 107 {
 108   //
 109   // Copy constructor
 110   //
 111 }
 112
 113 //_____________________________________________________________________________
 114 AliTRDtracker::AliTRDtracker(const TFile *geomfile)
 115   :AliTracker()
 116   ,fGeom(0)
 117   ,fNclusters(0)
 118   ,fClusters(new TObjArray(2000))
 119   ,fNseeds(0)
 120   ,fSeeds(new TObjArray(2000))
 121   ,fNtracks(0)
 122   ,fTracks(new TObjArray(1000))
 123   ,fTimeBinsPerPlane(0)
 124   ,fAddTRDseeds(kFALSE)
 125   ,fNoTilt(kFALSE)
 126   ,fDebugStreamer(0)
 127 {
 128   //
 129   //  Main constructor
 130   //
 131
 132   TDirectory *savedir = gDirectory;
 133   TFile      *in      = (TFile *) geomfile;
 134
 135   if (!in->IsOpen()) {
 136     AliWarning("geometry file is not open!\n");
 137     AliWarning("FULL TRD geometry and DEFAULT TRD parameter will be used\n");
 138   }
 139   else {
 140     in->cd();
 141     fGeom = (AliTRDgeometry *) in->Get("TRDgeometry");
 142   }
 143
 144   if (!fGeom) {
 145     AliWarning("Cannot find TRD geometry!\n");
 146     fGeom = new AliTRDgeometry();
 147   }
 148   fGeom->ReadGeoMatrices();
 149
 150   savedir->cd();
 151
 152   for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 153     Int_t trS   = CookSectorIndex(geomS);
 154     fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
 155     for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
 156       fHoles[icham][trS] = fGeom->IsHole(0,icham,geomS);
 157     }
 158   }
 159
 160   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
 161   Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
 162   if (tiltAngle < 0.1) {
 163     fNoTilt = kTRUE;
 164   }
 165
 166   fTimeBinsPerPlane =  AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
 167
 168   fDebugStreamer    = new TTreeSRedirector("TRDdebug.root");
 169
 170   savedir->cd();
 171
 172 }
 173
 174 //_____________________________________________________________________________
 175 AliTRDtracker::~AliTRDtracker()
 176 {
 177   //
 178   // Destructor of AliTRDtracker
 179   //
 180
 181   if (fClusters) {
 182     fClusters->Delete();
 183     delete fClusters;
 184   }
 185
 186   if (fTracks) {
 187     fTracks->Delete();
 188     delete fTracks;
 189   }
 190
 191   if (fSeeds) {
 192     fSeeds->Delete();
 193     delete fSeeds;
 194   }
 195
 196   delete fGeom;
 197
 198   for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 199     delete fTrSec[geomS];
 200   }
 201
 202   if (fDebugStreamer) {
 203     delete fDebugStreamer;
 204   }
 205
 206 }
 207
 208 //_____________________________________________________________________________
 209 Int_t  AliTRDtracker::LocalToGlobalID(Int_t lid)
 210 {
 211   //
 212   // Transform internal TRD ID to global detector ID
 213   //
 214
 215   Int_t  isector  = fGeom->GetSector(lid);
 216   Int_t  ichamber = fGeom->GetChamber(lid);
 217   Int_t  iplan    = fGeom->GetPlane(lid);
 218
 219   AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
 220   switch (iplan) {
 221   case 0:
 222     iLayer = AliAlignObj::kTRD1;
 223     break;
 224   case 1:
 225     iLayer = AliAlignObj::kTRD2;
 226     break;
 227   case 2:
 228     iLayer = AliAlignObj::kTRD3;
 229     break;
 230   case 3:
 231     iLayer = AliAlignObj::kTRD4;
 232     break;
 233   case 4:
 234     iLayer = AliAlignObj::kTRD5;
 235     break;
 236   case 5:
 237     iLayer = AliAlignObj::kTRD6;
 238     break;
 239   };
 240
 241   Int_t    modId = isector * fGeom->Ncham() + ichamber;
 242   UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
 243
 244   return volid;
 245
 246 }
 247
 248 //_____________________________________________________________________________
 249 Int_t  AliTRDtracker::GlobalToLocalID(Int_t gid)
 250 {
 251   //
 252   // Transform global detector ID to local detector ID
 253   //
 254
 255   Int_t modId    = 0;
 256   AliAlignObj::ELayerID layerId = AliAlignObj::VolUIDToLayer(gid,modId);
 257
 258   Int_t isector  = modId / fGeom->Ncham();
 259   Int_t ichamber = modId % fGeom->Ncham();
 260   Int_t iLayer   = -1;
 261
 262   switch (layerId) {
 263   case AliAlignObj::kTRD1:
 264     iLayer = 0;
 265     break;
 266   case AliAlignObj::kTRD2:
 267     iLayer = 1;
 268     break;
 269   case AliAlignObj::kTRD3:
 270     iLayer = 2;
 271     break;
 272   case AliAlignObj::kTRD4:
 273     iLayer = 3;
 274     break;
 275   case AliAlignObj::kTRD5:
 276     iLayer = 4;
 277     break;
 278   case AliAlignObj::kTRD6:
 279     iLayer = 5;
 280     break;
 281   default:
 282     iLayer =-1;
 283   }
 284
 285   if (iLayer < 0) {
 286     return -1;
 287   }
 288
 289   Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
 290
 291   return lid;
 292
 293 }
 294
 295 //_____________________________________________________________________________
 296 Bool_t  AliTRDtracker::Transform(AliTRDcluster *cluster)
 297 {
 298   //
 299   // Transform something ... whatever ...
 300   //
 301
 302   // Magic constants for geo manager transformation
 303   const Double_t kX0shift  = 2.52;
 304   const Double_t kX0shift5 = 3.05;
 305
 306   //
 307   // Apply alignment and calibration to transform cluster
 308   //
 309   Int_t detector = cluster->GetDetector();
 310   Int_t plane    = fGeom->GetPlane(cluster->GetDetector());
 311   Int_t chamber  = fGeom->GetChamber(cluster->GetDetector());
 312   Int_t sector   = fGeom->GetSector(cluster->GetDetector());
 313
 314   Double_t dxAmp  = (Double_t) fGeom->CamHght();               // Amplification region
 315   Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.0); // Drift distance
 316
 317   //
 318   // ExB correction
 319   //
 320   Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
 321   Double_t exB    = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
 322
 323   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
 324   AliTRDpadPlane    *padPlane    = commonParam->GetPadPlane(plane,chamber);
 325   Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
 326   Double_t localPos[3];
 327   Double_t localPosTracker[3];
 328   localPos[0] = -cluster->GetX();
 329   localPos[1] =  cluster->GetY() - driftX * exB;
 330   localPos[2] =  cluster->GetZ() - zshiftIdeal;
 331
 332   cluster->SetY(cluster->GetY() - driftX*exB);
 333   Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
 334   cluster->SetX(xplane- cluster->GetX());
 335
 336   TGeoHMatrix *matrix =  fGeom->GetCorrectionMatrix(cluster->GetDetector());
 337   if (!matrix) {
 338     // No matrix found - if somebody used geometry with holes
 339     AliError("Invalid Geometry - Default Geometry used\n");
 340     return kTRUE;
 341   }
 342   matrix->LocalToMaster(localPos,localPosTracker);
 343
 344   if (AliTRDReconstructor::StreamLevel() > 1) {
 345     (* fDebugStreamer) << "Transform"
 346                        << "Cl.="      << cluster
 347                        << "matrix.="  << matrix
 348                        << "Detector=" << detector
 349                        << "Sector="   << sector
 350                        << "Plane="    << plane
 351                        << "Chamber="  << chamber
 352                        << "lx0="      << localPosTracker[0]
 353                        << "ly0="      << localPosTracker[1]
 354                        << "lz0="      << localPosTracker[2]
 355                        << "\n";
 356   }
 357
 358   if (plane == 5) {
 359      cluster->SetX(localPosTracker[0]+kX0shift5);
 360   }
 361   else {
 362     cluster->SetX(localPosTracker[0]+kX0shift);
 363   }
 364
 365   cluster->SetY(localPosTracker[1]);
 366   cluster->SetZ(localPosTracker[2]);
 367
 368   return kTRUE;
 369
 370 }
 371
 372 //_____________________________________________________________________________
 373 // Bool_t  AliTRDtracker::Transform(AliTRDcluster *cluster)
 374 //{
 375 //   //
 376 //   // Is this still needed ????
 377 //   //
 378 //   const Double_t kDriftCorrection  = 1.01;                 // drift coeficient correction
 379 //   const Double_t kTime0Cor         = 0.32;                 // time0 correction
 380 //   //
 381 //   const Double_t kX0shift           = 2.52;
 382 //   const Double_t kX0shift5          = 3.05;
 383
 384 //   //
 385 //   // apply alignment and calibration to transform cluster
 386 //   //
 387 //   //
 388 //   Int_t detector = cluster->GetDetector();
 389 //   Int_t plane   = fGeom->GetPlane(cluster->GetDetector());
 390 //   Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
 391 //   Int_t sector  = fGeom->GetSector(cluster->GetDetector());
 392
 393 //   Double_t dxAmp  = (Double_t) fGeom->CamHght();          // Amplification region
 394 //   Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.);  // drift distance
 395 //   //
 396 //   // ExB correction
 397 //   //
 398 //   Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
 399 //   Double_t exB =   AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
 400 //   //
 401
 402 //   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
 403 //   AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
 404 //   Double_t zshiftIdeal  = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
 405 //   Double_t localPos[3], globalPos[3], localPosTracker[3], localPosTracker2[3];
 406 //   localPos[2] = -cluster->GetX();
 407 //   localPos[0] =  cluster->GetY() - driftX*exB;
 408 //   localPos[1] =  cluster->GetZ() -zshiftIdeal;
 409 //   TGeoHMatrix * matrix =  fGeom->GetGeoMatrix(cluster->GetDetector());
 410 //   matrix->LocalToMaster(localPos, globalPos);
 411
 412 //   Double_t sectorAngle = 20.*(sector%18)+10;
 413 //   TGeoHMatrix  rotSector;
 414 //   rotSector.RotateZ(sectorAngle);
 415 //   rotSector.LocalToMaster(globalPos, localPosTracker);
 416 //   //
 417 //   //
 418 //   TGeoHMatrix  matrix2(*matrix);
 419 //   matrix2.MultiplyLeft(&rotSector);
 420 //   matrix2.LocalToMaster(localPos,localPosTracker2);
 421 //   //
 422 //   //
 423 //   //
 424 //   cluster->SetY(cluster->GetY() - driftX*exB);
 425 //   Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
 426 //   cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor));
 427 //   (*fDebugStreamer)<<"Transform"<<
 428 //     "Cl.="<<cluster<<
 429 //     "matrix.="<<matrix<<
 430 //     "matrix2.="<<&matrix2<<
 431 //     "Detector="<<detector<<
 432 //     "Sector="<<sector<<
 433 //     "Plane="<<plane<<
 434 //     "Chamber="<<chamber<<
 435 //     "lx0="<<localPosTracker[0]<<
 436 //     "ly0="<<localPosTracker[1]<<
 437 //     "lz0="<<localPosTracker[2]<<
 438 //     "lx2="<<localPosTracker2[0]<<
 439 //     "ly2="<<localPosTracker2[1]<<
 440 //     "lz2="<<localPosTracker2[2]<<
 441 //     "\n";
 442 //   //
 443 //   if (plane==5)
 444 //      cluster->SetX(localPosTracker[0]+kX0shift5);
 445 //   else
 446 //     cluster->SetX(localPosTracker[0]+kX0shift);
 447
 448 //   cluster->SetY(localPosTracker[1]);
 449 //   cluster->SetZ(localPosTracker[2]);
 450 //   return kTRUE;
 451 // }
 452
 453 //_____________________________________________________________________________
 454 Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track)
 455 {
 456   //
 457   // Rotates the track when necessary
 458   //
 459
 460   Double_t alpha = AliTRDgeometry::GetAlpha();
 461   Double_t y     = track->GetY();
 462   Double_t ymax  = track->GetX()*TMath::Tan(0.5*alpha);
 463
 464   // Is this still needed ????
 465   //Int_t ns = AliTRDgeometry::kNsect;
 466   //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
 467
 468   if      (y >  ymax) {
 469     //s = (s+1) % ns;
 470     if (!track->Rotate( alpha)) {
 471       return kFALSE;
 472     }
 473   }
 474   else if (y < -ymax) {
 475     //s = (s-1+ns) % ns;
 476     if (!track->Rotate(-alpha)) {
 477       return kFALSE;
 478     }
 479   }
 480
 481   return kTRUE;
 482
 483 }
 484
 485 //_____________________________________________________________________________
 486 AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack *track, Int_t plane
 487                                        , Int_t timebin, UInt_t &index)
 488 {
 489   //
 490   // Try to find cluster in the backup list
 491   //
 492
 493   AliTRDcluster *cl =0;
 494   Int_t *indexes = track->GetBackupIndexes();
 495
 496   for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
 497     if (indexes[i] == 0) {
 498       break;
 499     }
 500     AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
 501     if (!cli) {
 502       break;
 503     }
 504     if (cli->GetLocalTimeBin() != timebin) {
 505       continue;
 506     }
 507     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 508     if (iplane == plane) {
 509       cl    = cli;
 510       index = indexes[i];
 511       break;
 512     }
 513   }
 514
 515   return cl;
 516
 517 }
 518
 519 //_____________________________________________________________________________
 520 Int_t  AliTRDtracker::GetLastPlane(AliTRDtrack *track)
 521 {
 522   //
 523   // Return last updated plane
 524   //
 525
 526   Int_t  lastplane = 0;
 527   Int_t *indexes   = track->GetBackupIndexes();
 528
 529   for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
 530     AliTRDcluster *cli = (AliTRDcluster *) fClusters->UncheckedAt(indexes[i]);
 531     if (!cli) {
 532       break;
 533     }
 534     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 535     if (iplane > lastplane) {
 536       lastplane = iplane;
 537     }
 538   }
 539
 540   return lastplane;
 541
 542 }
 543
 544 //_____________________________________________________________________________
 545 Int_t AliTRDtracker::Clusters2Tracks(AliESD *event)
 546 {
 547   //
 548   // Finds tracks within the TRD. The ESD event is expected to contain seeds
 549   // at the outer part of the TRD. The seeds
 550   // are found within the TRD if fAddTRDseeds is TRUE.
 551   // The tracks are propagated to the innermost time bin
 552   // of the TRD and the ESD event is updated
 553   //
 554
 555   Int_t   timeBins = fTrSec[0]->GetNumberOfTimeBins();
 556   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 557   Int_t   nseed    = 0;
 558   Int_t   found    = 0;
 559   //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 560
 561   Int_t n = event->GetNumberOfTracks();
 562   for (Int_t i = 0; i < n; i++) {
 563
 564     AliESDtrack *seed = event->GetTrack(i);
 565     ULong_t status = seed->GetStatus();
 566     if ((status & AliESDtrack::kTRDout) == 0) {
 567       continue;
 568     }
 569     if ((status & AliESDtrack::kTRDin)  != 0) {
 570       continue;
 571     }
 572     nseed++;
 573
 574     AliTRDtrack *seed2 = new AliTRDtrack(*seed);
 575     //seed2->ResetCovariance();
 576     AliTRDtrack *pt    = new AliTRDtrack(*seed2,seed2->GetAlpha());
 577     AliTRDtrack &t     = *pt;
 578     FollowProlongation(t);
 579     if (t.GetNumberOfClusters() >= foundMin) {
 580       UseClusters(&t);
 581       CookLabel(pt,1 - fgkLabelFraction);
 582       //t.CookdEdx();
 583     }
 584     found++;
 585
 586     Double_t xTPC = 250.0;
 587     if (PropagateToX(t,xTPC,fgkMaxStep)) {
 588       seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
 589     }
 590     delete seed2;
 591     delete pt;
 592
 593   }
 594
 595   AliInfo(Form("Number of loaded seeds: %d",nseed));
 596   AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
 597   AliInfo(Form("Total number of found tracks: %d",found));
 598
 599   return 0;
 600
 601 }
 602
 603 //_____________________________________________________________________________
 604 Int_t AliTRDtracker::PropagateBack(AliESD *event)
 605 {
 606   //
 607   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
 608   // backpropagated by the TPC tracker. Each seed is first propagated
 609   // to the TRD, and then its prolongation is searched in the TRD.
 610   // If sufficiently long continuation of the track is found in the TRD
 611   // the track is updated, otherwise it's stored as originaly defined
 612   // by the TPC tracker.
 613   //
 614
 615   Int_t   found    = 0;
 616   Float_t foundMin = 20.0;
 617   Int_t   n        = event->GetNumberOfTracks();
 618
 619   // Sort tracks
 620   Float_t *quality = new Float_t[n];
 621   Int_t   *index   = new Int_t[n];
 622   for (Int_t i = 0; i < n; i++) {
 623     AliESDtrack *seed = event->GetTrack(i);
 624     Double_t covariance[15];
 625     seed->GetExternalCovariance(covariance);
 626     quality[i] = covariance[0]+covariance[2];
 627   }
 628   TMath::Sort(n,quality,index,kFALSE);
 629
 630   for (Int_t i = 0; i < n; i++) {
 631
 632     //AliESDtrack *seed = event->GetTrack(i);
 633     AliESDtrack *seed = event->GetTrack(index[i]);
 634
 635     ULong_t status = seed->GetStatus();
 636     if ((status & AliESDtrack::kTPCout) == 0) {
 637       continue;
 638     }
 639     if ((status & AliESDtrack::kTRDout) != 0) {
 640       continue;
 641     }
 642
 643     Int_t   lbl = seed->GetLabel();
 644     AliTRDtrack *track = new AliTRDtrack(*seed);
 645     track->SetSeedLabel(lbl);
 646     seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
 647     fNseeds++;
 648     Float_t p4  = track->GetC();
 649
 650     Int_t expectedClr = FollowBackProlongation(*track);
 651     if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
 652         (TMath::Abs(track->GetPt())                      > 0.8)) {
 653
 654       //
 655       // Make backup for back propagation
 656       //
 657
 658       Int_t foundClr = track->GetNumberOfClusters();
 659       if (foundClr >= foundMin) {
 660         track->CookdEdx();
 661         CookdEdxTimBin(*track);
 662         CookLabel(track,1 - fgkLabelFraction);
 663         if (track->GetBackupTrack()) {
 664           UseClusters(track->GetBackupTrack());
 665         }
 666
 667         // Sign only gold tracks
 668         if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
 669           if ((seed->GetKinkIndex(0)      ==   0) &&
 670               (TMath::Abs(track->GetPt()) <  1.5)) {
 671             UseClusters(track);
 672           }
 673         }
 674         Bool_t isGold = kFALSE;
 675
 676         // Full gold track
 677         if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
 678           //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 679           if (track->GetBackupTrack()) {
 680             seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 681           }
 682           isGold = kTRUE;
 683         }
 684
 685         // Almost gold track
 686         if ((!isGold)                                              &&
 687             (track->GetNCross()                              == 0) &&
 688             (track->GetChi2() / track->GetNumberOfClusters()  < 7)) {
 689           //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 690           if (track->GetBackupTrack()) {
 691             seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 692           }
 693           isGold = kTRUE;
 694         }
 695
 696         if ((!isGold) &&
 697             (track->GetBackupTrack())) {
 698           if ((track->GetBackupTrack()->GetNumberOfClusters()                                          > foundMin) &&
 699               ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) <        7)) {
 700             seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 701             isGold = kTRUE;
 702           }
 703         }
 704
 705         if ((track->StatusForTOF()                          >   0) &&
 706             (track->fNCross                                ==   0) &&
 707             (Float_t(track->fN)/Float_t(track->fNExpected)  > 0.4)) {
 708           //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 709         }
 710
 711       }
 712
 713     }
 714
 715     // Debug part of tracking
 716     TTreeSRedirector &cstream = *fDebugStreamer;
 717     Int_t eventNr = event->GetEventNumber();
 718     if (AliTRDReconstructor::StreamLevel() > 0) {
 719       if (track->GetBackupTrack()) {
 720         cstream << "Tracks"
 721                 << "EventNr="  << eventNr
 722                 << "ESD.="     << seed
 723                 << "trd.="     << track
 724                 << "trdback.=" << track->GetBackupTrack()
 725                 << "\n";
 726       }
 727       else {
 728         cstream << "Tracks"
 729                 << "EventNr="  << eventNr
 730                 << "ESD.="     << seed
 731                 << "trd.="     << track
 732                 << "trdback.=" << track
 733                 << "\n";
 734       }
 735     }
 736
 737     // Propagation to the TOF (I.Belikov)
 738     if (track->GetStop() == kFALSE) {
 739
 740       Double_t xtof  = 371.0;
 741       Double_t c2    = track->GetSnp() + track->GetC() * (xtof - track->GetX());
 742       if (TMath::Abs(c2) >= 0.99) {
 743         delete track;
 744         continue;
 745       }
 746       Double_t xTOF0 = 370.0;
 747       PropagateToX(*track,xTOF0,fgkMaxStep);
 748
 749       // Energy losses taken to the account - check one more time
 750       c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
 751       if (TMath::Abs(c2) >= 0.99) {
 752         delete track;
 753         continue;
 754       }
 755
 756       Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
 757       Double_t y;
 758       track->GetYAt(xtof,GetBz(),y);
 759       if      (y >  ymax) {
 760         if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
 761           delete track;
 762           continue;
 763         }
 764       }
 765       else if (y < -ymax) {
 766         if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
 767           delete track;
 768           continue;
 769         }
 770       }
 771
 772       if (track->PropagateTo(xtof)) {
 773         seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
 774         for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 775           for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
 776             seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
 777           }
 778           seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 779         }
 780         //seed->SetTRDtrack(new AliTRDtrack(*track));
 781         if (track->GetNumberOfClusters() > foundMin) {
 782           found++;
 783         }
 784       }
 785
 786     }
 787     else {
 788
 789       if ((track->GetNumberOfClusters() >              15) &&
 790           (track->GetNumberOfClusters() > 0.5*expectedClr)) {
 791         seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
 792         //seed->SetStatus(AliESDtrack::kTRDStop);
 793         for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 794           for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
 795             seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
 796           }
 797           seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 798         }
 799         //seed->SetTRDtrack(new AliTRDtrack(*track));
 800         found++;
 801       }
 802
 803     }
 804
 805     seed->SetTRDQuality(track->StatusForTOF());
 806     seed->SetTRDBudget(track->fBudget[0]);
 807
 808     delete track;
 809
 810   }
 811
 812   AliInfo(Form("Number of seeds: %d",fNseeds));
 813   AliInfo(Form("Number of back propagated TRD tracks: %d",found));
 814
 815   // New seeding
 816   if (AliTRDReconstructor::SeedingOn()) {
 817     MakeSeedsMI(3,5,event);
 818   }
 819
 820   fSeeds->Clear();
 821   fNseeds = 0;
 822
 823   delete [] index;
 824   delete [] quality;
 825
 826   return 0;
 827
 828 }
 829
 830 //_____________________________________________________________________________
 831 Int_t AliTRDtracker::RefitInward(AliESD *event)
 832 {
 833   //
 834   // Refits tracks within the TRD. The ESD event is expected to contain seeds
 835   // at the outer part of the TRD.
 836   // The tracks are propagated to the innermost time bin
 837   // of the TRD and the ESD event is updated
 838   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
 839   //
 840
 841   Int_t   timeBins = fTrSec[0]->GetNumberOfTimeBins();
 842   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 843   Int_t   nseed    = 0;
 844   Int_t   found    = 0;
 845   //Int_t innerTB    = fTrSec[0]->GetInnerTimeBin();
 846   AliTRDtrack seed2;
 847
 848   Int_t n = event->GetNumberOfTracks();
 849   for (Int_t i = 0; i < n; i++) {
 850
 851     AliESDtrack *seed = event->GetTrack(i);
 852     new (&seed2) AliTRDtrack(*seed);
 853     if (seed2.GetX() < 270.0) {
 854       seed->UpdateTrackParams(&seed2,AliESDtrack::kTRDbackup); // Backup TPC track - only update
 855       continue;
 856     }
 857
 858     ULong_t status = seed->GetStatus();
 859     if ((status & AliESDtrack::kTRDout) == 0) {
 860       continue;
 861     }
 862     if ((status & AliESDtrack::kTRDin)  != 0) {
 863       continue;
 864     }
 865     nseed++;
 866
 867     seed2.ResetCovariance(50.0);
 868
 869     AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
 870     Int_t *indexes2 = seed2.GetIndexes();
 871     for (Int_t i = 0; i < AliESDtrack::kNPlane;i++) {
 872       for (Int_t j = 0; j < AliESDtrack::kNSlice;j++) {
 873         pt->SetPIDsignals(seed2.GetPIDsignals(i,j),i,j);
 874       }
 875       pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
 876     }
 877
 878     Int_t *indexes3 = pt->GetBackupIndexes();
 879     for (Int_t i = 0; i < 200;i++) {
 880       if (indexes2[i] == 0) {
 881         break;
 882       }
 883       indexes3[i] = indexes2[i];
 884     }
 885
 886     //AliTRDtrack *pt = seed2;
 887     AliTRDtrack &t = *pt;
 888     FollowProlongation(t);
 889     if (t.GetNumberOfClusters() >= foundMin) {
 890       //UseClusters(&t);
 891       //CookLabel(pt, 1-fgkLabelFraction);
 892       t.CookdEdx();
 893       CookdEdxTimBin(t);
 894     }
 895     found++;
 896
 897     Double_t xTPC = 250.0;
 898     if (PropagateToX(t,xTPC,fgkMaxStep)) {
 899       seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
 900       for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 901         for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
 902           seed->SetTRDsignals(pt->GetPIDsignals(i,j),i,j);
 903         }
 904         seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
 905       }
 906     }
 907     else {
 908       // If not prolongation to TPC - propagate without update
 909       AliTRDtrack *seed2 = new AliTRDtrack(*seed);
 910       seed2->ResetCovariance(5.0);
 911       AliTRDtrack *pt2   = new AliTRDtrack(*seed2,seed2->GetAlpha());
 912       delete seed2;
 913       if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
 914         pt2->CookdEdx( );
 915         CookdEdxTimBin(*pt2);
 916         seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
 917         for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 918           for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
 919             seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
 920           }
 921           seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
 922         }
 923       }
 924       delete pt2;
 925     }
 926
 927     delete pt;
 928
 929   }
 930
 931   AliInfo(Form("Number of loaded seeds: %d",nseed));
 932   AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
 933
 934   return 0;
 935
 936 }
 937
 938 //_____________________________________________________________________________
 939 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
 940 {
 941   //
 942   // Starting from current position on track=t this function tries
 943   // to extrapolate the track up to timeBin=0 and to confirm prolongation
 944   // if a close cluster is found. Returns the number of clusters
 945   // expected to be found in sensitive layers
 946   // GeoManager used to estimate mean density
 947   //
 948
 949   Int_t    sector;
 950   Int_t    lastplane = GetLastPlane(&t);
 951   Double_t radLength = 0.0;
 952   Double_t rho       = 0.0;
 953   Int_t    expectedNumberOfClusters = 0;
 954
 955   for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
 956
 957     Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
 958     Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
 959
 960     //
 961     // Propagate track close to the plane if neccessary
 962     //
 963     Double_t currentx  = fTrSec[0]->GetLayer(rowlast)->GetX();
 964     if (currentx < (-fgkMaxStep + t.GetX())) {
 965       // Propagate closer to chamber - safety space fgkMaxStep
 966       if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
 967         break;
 968       }
 969     }
 970
 971     if (!AdjustSector(&t)) {
 972       break;
 973     }
 974
 975     //
 976     // Get material budget
 977     //
 978     Double_t xyz0[3];
 979     Double_t xyz1[3];
 980     Double_t param[7];
 981     Double_t x;
 982     Double_t y;
 983     Double_t z;
 984
 985     // Starting global position
 986     t.GetXYZ(xyz0);
 987     // End global position
 988     x = fTrSec[0]->GetLayer(row0)->GetX();
 989     if (!t.GetProlongation(x,y,z)) {
 990       break;
 991     }
 992     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
 993     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
 994     xyz1[2] =  z;
 995     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 996     rho       = param[0];
 997     radLength = param[1]; // Get mean propagation parameters
 998
 999     //
1000     // Propagate and update
1001     //
1002     sector = t.GetSector();
1003     //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
1004     for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
1005
1006       Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
1007       expectedNumberOfClusters++;
1008       t.fNExpected++;
1009       if (t.GetX() > 345.0) {
1010         t.fNExpectedLast++;
1011       }
1012       AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
1013       AliTRDcluster          *cl      = 0;
1014       UInt_t   index   = 0;
1015       Double_t maxChi2 = fgkMaxChi2;
1016       x = timeBin.GetX();
1017
1018       if (timeBin) {
1019
1020         AliTRDcluster *cl0 = timeBin[0];
1021         if (!cl0) {
1022           // No clusters in given time bin
1023           continue;
1024         }
1025
1026         Int_t plane   = fGeom->GetPlane(cl0->GetDetector());
1027         if (plane > lastplane) {
1028           continue;
1029         }
1030
1031         Int_t timebin = cl0->GetLocalTimeBin();
1032         AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index);
1033
1034         if (cl2) {
1035
1036           cl = cl2;
1037           //Double_t h01 = GetTiltFactor(cl);    //I.B's fix
1038           //maxChi2=t.GetPredictedChi2(cl,h01);
1039
1040         }
1041         if (cl) {
1042
1043           //if (cl->GetNPads()<5)
1044           Double_t dxsample = timeBin.GetdX();
1045           t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1046           Double_t h01      = GetTiltFactor(cl);
1047           Int_t    det      = cl->GetDetector();
1048           Int_t    plane    = fGeom->GetPlane(det);
1049           if (t.GetX() > 345.0) {
1050             t.fNLast++;
1051             t.fChi2Last += maxChi2;
1052           }
1053
1054           Double_t xcluster = cl->GetX();
1055           t.PropagateTo(xcluster,radLength,rho);
1056
1057           if (!AdjustSector(&t)) {
1058             break;           //I.B's fix
1059           }
1060           maxChi2 = t.GetPredictedChi2(cl,h01);
1061
1062           if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1063             // ????
1064           }
1065
1066         }
1067
1068       }
1069
1070     }
1071
1072   }
1073
1074   return expectedNumberOfClusters;
1075
1076 }
1077
1078 //_____________________________________________________________________________
1079 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
1080 {
1081   //
1082   // Starting from current radial position of track <t> this function
1083   // extrapolates the track up to outer timebin and in the sensitive
1084   // layers confirms prolongation if a close cluster is found.
1085   // Returns the number of clusters expected to be found in sensitive layers
1086   // Use GEO manager for material Description
1087   //
1088
1089   Int_t    sector;
1090   Int_t    clusters[1000];
1091   Double_t radLength = 0.0;
1092   Double_t rho       = 0.0;
1093   Int_t    expectedNumberOfClusters = 0;
1094   Float_t  ratio0    = 0.0;
1095   AliTRDtracklet tracklet;
1096
1097   // Calibration fill 2D
1098   AliTRDCalibra *calibra = AliTRDCalibra::Instance();
1099   if (!calibra) {
1100     AliInfo("Could not get Calibra instance\n");
1101   }
1102   if (calibra->GetMItracking()) {
1103     calibra->Resettrack();
1104   }
1105
1106   for (Int_t i = 0; i < 1000; i++) {
1107     clusters[i] = -1;
1108   }
1109
1110   for (Int_t iplane = 0; iplane < AliESDtrack::kNPlane; iplane++) {
1111
1112     Int_t    row0     = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
1113     Int_t    rowlast  = GetGlobalTimeBin(0,iplane,0);
1114     Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
1115     if (currentx < t.GetX()) {
1116       continue;
1117     }
1118
1119     //
1120     // Propagate closer to chamber if neccessary
1121     //
1122     if (currentx > (fgkMaxStep + t.GetX())) {
1123       if (!PropagateToX(t,currentx-fgkMaxStep,fgkMaxStep)) {
1124         break;
1125       }
1126     }
1127     if (!AdjustSector(&t)) {
1128       break;
1129     }
1130     if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
1131       break;
1132     }
1133
1134     //
1135     // Get material budget inside of chamber
1136     //
1137     Double_t xyz0[3];
1138     Double_t xyz1[3];
1139     Double_t param[7];
1140     Double_t x;
1141     Double_t y;
1142     Double_t z;
1143     // Starting global position
1144     t.GetXYZ(xyz0);
1145     // End global position
1146     x = fTrSec[0]->GetLayer(rowlast)->GetX();
1147     if (!t.GetProlongation(x,y,z)) {
1148       break;
1149     }
1150     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1151     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1152     xyz1[2] =  z;
1153     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1154     rho       = param[0];
1155     radLength = param[1]; // Get mean propagation parameters
1156
1157     //
1158     // Find clusters
1159     //
1160     sector = t.GetSector();
1161     Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
1162     if (tracklet.GetN() < GetTimeBinsPerPlane()/3) {
1163       continue;
1164     }
1165
1166     //
1167     // Propagate and update track
1168     //
1169     for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
1170
1171       Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
1172       expectedNumberOfClusters++;
1173       t.fNExpected++;
1174       if (t.GetX() > 345.0) {
1175         t.fNExpectedLast++;
1176       }
1177       AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
1178       AliTRDcluster          *cl      = 0;
1179       UInt_t   index   = 0;
1180       Double_t maxChi2 = fgkMaxChi2;
1181       x = timeBin.GetX();
1182
1183       if (timeBin) {
1184
1185         if (clusters[ilayer] > 0) {
1186           index = clusters[ilayer];
1187           cl    = (AliTRDcluster *)GetCluster(index);
1188           //Double_t h01 = GetTiltFactor(cl);   // I.B's fix
1189           //maxChi2=t.GetPredictedChi2(cl,h01); //
1190         }
1191
1192         if (cl) {
1193
1194           //if (cl->GetNPads() < 5)
1195           Double_t dxsample = timeBin.GetdX();
1196           t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1197           Double_t h01      = GetTiltFactor(cl);
1198           Int_t    det      = cl->GetDetector();
1199           Int_t    plane    = fGeom->GetPlane(det);
1200           if (t.GetX() > 345.0) {
1201             t.fNLast++;
1202             t.fChi2Last += maxChi2;
1203           }
1204           Double_t xcluster = cl->GetX();
1205           t.PropagateTo(xcluster,radLength,rho);
1206           maxChi2 = t.GetPredictedChi2(cl,h01);
1207
1208           if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1209             if (!t.Update(cl,maxChi2,index,h01)) {
1210               // ????
1211             }
1212           }
1213
1214           if (calibra->GetMItracking()) {
1215             calibra->UpdateHistograms(cl,&t);
1216           }
1217
1218           // Reset material budget if 2 consecutive gold
1219           if (plane > 0) {
1220             if (t.fTracklets[plane].GetN() + t.fTracklets[plane-1].GetN() > 20) {
1221               t.fBudget[2] = 0;
1222             }
1223           }
1224
1225         }
1226
1227       }
1228
1229     }
1230
1231     ratio0 = ncl / Float_t(fTimeBinsPerPlane);
1232     Float_t ratio1 = Float_t(t.fN+1) / Float_t(t.fNExpected+1.0);
1233     if ((tracklet.GetChi2()     <  18.0) &&
1234         (ratio0                 >   0.8) &&
1235         (ratio1                 >   0.6) &&
1236         (ratio0+ratio1          >   1.5) &&
1237         (t.GetNCross()          ==    0) &&
1238         (TMath::Abs(t.GetSnp()) <  0.85) &&
1239         (t.fN                   >    20)){
1240       t.MakeBackupTrack(); // Make backup of the track until is gold
1241     }
1242
1243   }
1244
1245   return expectedNumberOfClusters;
1246
1247 }
1248
1249 //_____________________________________________________________________________
1250 Int_t AliTRDtracker::PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep)
1251 {
1252   //
1253   // Starting from current radial position of track <t> this function
1254   // extrapolates the track up to radial position <xToGo>.
1255   // Returns 1 if track reaches the plane, and 0 otherwise
1256   //
1257
1258   const Double_t kEpsilon = 0.00001;
1259   //Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
1260   Double_t xpos = t.GetX();
1261   Double_t dir  = (xpos<xToGo) ? 1.0 : -1.0;
1262
1263   while (((xToGo-xpos)*dir) > kEpsilon) {
1264
1265     Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1266
1267     Double_t xyz0[3];
1268     Double_t xyz1[3];
1269     Double_t param[7];
1270     Double_t x;
1271     Double_t y;
1272     Double_t z;
1273     // Starting global position
1274     t.GetXYZ(xyz0);
1275     x = xpos + step;
1276
1277     if (!t.GetProlongation(x,y,z)) {
1278       return 0; // No prolongation
1279     }
1280
1281     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1282     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1283     xyz1[2] =  z;
1284
1285     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1286     if (!t.PropagateTo(x,param[1],param[0])) {
1287       return 0;
1288     }
1289     AdjustSector(&t);
1290     xpos = t.GetX();
1291
1292   }
1293
1294   return 1;
1295
1296 }
1297
1298 //_____________________________________________________________________________
1299 Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1300 {
1301   //
1302   // Fills clusters into TRD tracking_sectors
1303   // Note that the numbering scheme for the TRD tracking_sectors
1304   // differs from that of TRD sectors
1305   //
1306
1307   if (ReadClusters(fClusters,cTree)) {
1308      AliError("Problem with reading the clusters !");
1309      return 1;
1310   }
1311   Int_t ncl = fClusters->GetEntriesFast();
1312   fNclusters = ncl;
1313   AliInfo(Form("Sorting %d clusters",ncl));
1314
1315   UInt_t index;
1316   for (Int_t ichamber = 0; ichamber <  5; ichamber++) {
1317     for (Int_t isector  = 0; isector  < 18;  isector++) {
1318       fHoles[ichamber][isector] = kTRUE;
1319     }
1320   }
1321
1322   while (ncl--) {
1323
1324     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(ncl);
1325     Int_t detector       = c->GetDetector();
1326     Int_t localTimeBin   = c->GetLocalTimeBin();
1327     Int_t sector         = fGeom->GetSector(detector);
1328     Int_t plane          = fGeom->GetPlane(detector);
1329     Int_t trackingSector = CookSectorIndex(sector);
1330
1331     if (c->GetLabel(0) > 0) {
1332       Int_t chamber = fGeom->GetChamber(detector);
1333       fHoles[chamber][trackingSector] = kFALSE;
1334     }
1335
1336     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1337     if (gtb < 0) {
1338       continue;
1339     }
1340     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1341
1342     index = ncl;
1343
1344     // Apply pos correction
1345     Transform(c);
1346     fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
1347
1348   }
1349
1350   return 0;
1351
1352 }
1353
1354 //_____________________________________________________________________________
1355 void AliTRDtracker::UnloadClusters()
1356 {
1357   //
1358   // Clears the arrays of clusters and tracks. Resets sectors and timebins
1359   //
1360
1361   Int_t i;
1362   Int_t nentr;
1363
1364   nentr = fClusters->GetEntriesFast();
1365   for (i = 0; i < nentr; i++) {
1366     delete fClusters->RemoveAt(i);
1367   }
1368   fNclusters = 0;
1369
1370   nentr = fSeeds->GetEntriesFast();
1371   for (i = 0; i < nentr; i++) {
1372     delete fSeeds->RemoveAt(i);
1373   }
1374
1375   nentr = fTracks->GetEntriesFast();
1376   for (i = 0; i < nentr; i++) {
1377     delete fTracks->RemoveAt(i);
1378   }
1379
1380   Int_t nsec = AliTRDgeometry::kNsect;
1381   for (i = 0; i < nsec; i++) {
1382     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1383       fTrSec[i]->GetLayer(pl)->Clear();
1384     }
1385   }
1386
1387 }
1388
1389 //_____________________________________________________________________________
1390 void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd)
1391 {
1392   //
1393   // Creates  seeds using clusters between  position inner plane  and outer plane
1394   //
1395
1396   const Double_t kMaxTheta =  1.0;
1397   const Double_t kMaxPhi   =  2.0;
1398
1399   const Double_t kRoad0y   =  6.0; // Road for middle cluster
1400   const Double_t kRoad0z   =  8.5; // Road for middle cluster
1401
1402   const Double_t kRoad1y   =  2.0; // Road in y for seeded cluster
1403   const Double_t kRoad1z   = 20.0; // Road in z for seeded cluster
1404
1405   const Double_t kRoad2y   =  3.0; // Road in y for extrapolated cluster
1406   const Double_t kRoad2z   = 20.0; // Road in z for extrapolated cluster
1407   const Int_t    kMaxSeed  = 3000;
1408
1409   Int_t maxSec = AliTRDgeometry::kNsect;
1410
1411   // Linear fitters in planes
1412   TLinearFitter fitterTC(2,"hyp2");  // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1413   TLinearFitter fitterT2(4,"hyp4");  // Fitting with tilting pads - kz not fixed
1414   fitterTC.StoreData(kTRUE);
1415   fitterT2.StoreData(kTRUE);
1416   AliRieman rieman(1000);            // Rieman fitter
1417   AliRieman rieman2(1000);           // Rieman fitter
1418
1419   // Find the maximal and minimal layer for the planes
1420   Int_t layers[6][2];
1421   AliTRDpropagationLayer *reflayers[6];
1422   for (Int_t i = 0; i < 6; i++) {
1423     layers[i][0] = 10000;
1424     layers[i][1] = 0;
1425   }
1426   for (Int_t ns = 0; ns < maxSec; ns++) {
1427     for (Int_t ilayer = 0; ilayer < fTrSec[ns]->GetNumberOfLayers(); ilayer++) {
1428       AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(ilayer));
1429       if (layer == 0) {
1430         continue;
1431       }
1432       Int_t det   = layer[0]->GetDetector();
1433       Int_t plane = fGeom->GetPlane(det);
1434       if (ilayer < layers[plane][0]) {
1435         layers[plane][0] = ilayer;
1436       }
1437       if (ilayer > layers[plane][1]) {
1438         layers[plane][1] = ilayer;
1439       }
1440     }
1441   }
1442
1443   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
1444   Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
1445   Double_t hL[6];                                                // Tilting angle
1446   Double_t xcl[6];                                               // X - position of reference cluster
1447   Double_t ycl[6];                                               // Y - position of reference cluster
1448   Double_t zcl[6];                                               // Z - position of reference cluster
1449
1450   AliTRDcluster *cl[6] = {    0,    0,    0,    0,    0,    0 }; // Seeding clusters
1451   Float_t padlength[6] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 }; // Current pad-length
1452
1453   Double_t chi2R  = 0.0;
1454   Double_t chi2Z  = 0.0;
1455   Double_t chi2RF = 0.0;
1456   Double_t chi2ZF = 0.0;
1457
1458   Int_t nclusters;                                               // Total number of clusters
1459   for (Int_t i = 0; i < 6; i++) {
1460     hL[i] = h01;
1461     if (i%2==1) {
1462       hL[i]*=-1.0;
1463     }
1464   }
1465
1466   // Registered seed
1467   AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1468   AliTRDseed *seed[kMaxSeed];
1469   for (Int_t iseed = 0; iseed < kMaxSeed; iseed++) {
1470     seed[iseed]= &pseed[iseed*6];
1471   }
1472   AliTRDseed *cseed = seed[0];
1473
1474   Double_t seedquality[kMaxSeed];
1475   Double_t seedquality2[kMaxSeed];
1476   Double_t seedparams[kMaxSeed][7];
1477   Int_t    seedlayer[kMaxSeed];
1478   Int_t    registered = 0;
1479   Int_t    sort[kMaxSeed];
1480
1481   //
1482   // Seeding part
1483   //
1484   for (Int_t ns = 0; ns < maxSec; ns++) { // Loop over sectors
1485   //for (Int_t ns = 0; ns < 5; ns++) { // Loop over sectors
1486
1487     registered   = 0; // Reset registerd seed counter
1488     cseed        = seed[registered];
1489     Float_t iter = 0.0;
1490
1491     for (Int_t sLayer = 2; sLayer >= 0; sLayer--) {
1492       //for (Int_t dseed = 5; dseed < 15; dseed += 3) {
1493
1494       iter += 1.0;
1495       Int_t dseed = 5 + Int_t(iter) * 3;
1496
1497       // Initialize seeding layers
1498       for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
1499         reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1500         xcl[ilayer]       = reflayers[ilayer]->GetX();
1501       }
1502
1503       Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2]) * 0.5;
1504       AliTRDpropagationLayer &layer0 = *reflayers[sLayer+0];
1505       AliTRDpropagationLayer &layer1 = *reflayers[sLayer+1];
1506       AliTRDpropagationLayer &layer2 = *reflayers[sLayer+2];
1507       AliTRDpropagationLayer &layer3 = *reflayers[sLayer+3];
1508
1509       Int_t maxn3 = layer3;
1510       for (Int_t icl3 = 0; icl3 < maxn3; icl3++) {
1511
1512         AliTRDcluster *cl3 = layer3[icl3];
1513         if (!cl3) {
1514           continue;
1515         }
1516         padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2() * 12.0);
1517         ycl[sLayer+3] = cl3->GetY();
1518         zcl[sLayer+3] = cl3->GetZ();
1519         Float_t yymin0 = ycl[sLayer+3] - 1.0 - kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1520         Float_t yymax0 = ycl[sLayer+3] + 1.0 + kMaxPhi * (xcl[sLayer+3]-xcl[sLayer+0]);
1521         Int_t   maxn0  = layer0;
1522
1523         for (Int_t icl0 = layer0.Find(yymin0); icl0 < maxn0; icl0++) {
1524
1525           AliTRDcluster *cl0 = layer0[icl0];
1526           if (!cl0) {
1527             continue;
1528           }
1529           if (cl3->IsUsed() && cl0->IsUsed()) {
1530             continue;
1531           }
1532           ycl[sLayer+0] = cl0->GetY();
1533           zcl[sLayer+0] = cl0->GetZ();
1534           if (ycl[sLayer+0] > yymax0) {
1535             break;
1536           }
1537           Double_t tanphi     = (ycl[sLayer+3]-ycl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1538           if (TMath::Abs(tanphi)   > kMaxPhi) {
1539             continue;
1540           }
1541           Double_t tantheta   = (zcl[sLayer+3]-zcl[sLayer+0]) / (xcl[sLayer+3]-xcl[sLayer+0]);
1542           if (TMath::Abs(tantheta) > kMaxTheta) {
1543             continue;
1544           }
1545           padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2() * 12.0);
1546
1547           // Expected position in 1 layer
1548           Double_t y1exp = ycl[sLayer+0] + (tanphi)   * (xcl[sLayer+1]-xcl[sLayer+0]);
1549           Double_t z1exp = zcl[sLayer+0] + (tantheta) * (xcl[sLayer+1]-xcl[sLayer+0]);
1550           Float_t yymin1 = y1exp - kRoad0y - tanphi;
1551           Float_t yymax1 = y1exp + kRoad0y + tanphi;
1552           Int_t   maxn1  = layer1;
1553
1554           for (Int_t icl1 = layer1.Find(yymin1); icl1 < maxn1; icl1++) {
1555
1556             AliTRDcluster *cl1 = layer1[icl1];
1557             if (!cl1) {
1558               continue;
1559             }
1560             Int_t nusedCl = 0;
1561             if (cl3->IsUsed()) nusedCl++;
1562             if (cl0->IsUsed()) nusedCl++;
1563             if (cl1->IsUsed()) nusedCl++;
1564             if (nusedCl > 1) {
1565               continue;
1566             }
1567             ycl[sLayer+1] = cl1->GetY();
1568             zcl[sLayer+1] = cl1->GetZ();
1569             if (ycl[sLayer+1] > yymax1) {
1570               break;
1571             }
1572             if (TMath::Abs(ycl[sLayer+1]-y1exp) > kRoad0y+tanphi) {
1573               continue;
1574             }
1575             if (TMath::Abs(zcl[sLayer+1]-z1exp) > kRoad0z) {
1576               continue;
1577             }
1578             padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2() * 12.0);
1579
1580             Double_t y2exp  = ycl[sLayer+0]+(tanphi)   * (xcl[sLayer+2]-xcl[sLayer+0]) + (ycl[sLayer+1]-y1exp);
1581             Double_t z2exp  = zcl[sLayer+0]+(tantheta) * (xcl[sLayer+2]-xcl[sLayer+0]);
1582             Int_t    index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,kRoad1z);
1583             if (index2 <= 0) {
1584               continue;
1585             }
1586             AliTRDcluster *cl2 = (AliTRDcluster *) GetCluster(index2);
1587             padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2() * 12.0);
1588             ycl[sLayer+2]       = cl2->GetY();
1589             zcl[sLayer+2]       = cl2->GetZ();
1590             if (TMath::Abs(cl2->GetZ()-z2exp) > kRoad0z) {
1591               continue;
1592             }
1593
1594             rieman.Reset();
1595             rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1596             rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1597             rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1598             rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1599             rieman.Update();
1600
1601             // Reset fitter
1602             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1603               cseed[iLayer].Reset();
1604             }
1605             chi2Z = 0.0;
1606             chi2R = 0.0;
1607
1608             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1609               cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1610               chi2Z += (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer])
1611                      * (cseed[sLayer+iLayer].GetZref(0)- zcl[sLayer+iLayer]);
1612               cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1613               cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1614               chi2R += (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer])
1615                      * (cseed[sLayer+iLayer].GetYref(0)- ycl[sLayer+iLayer]);
1616               cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1617             }
1618             if (TMath::Sqrt(chi2R) > 1.0/iter) {
1619               continue;
1620             }
1621             if (TMath::Sqrt(chi2Z) > 7.0/iter) {
1622               continue;
1623             }
1624
1625             Float_t minmax[2] = { -100.0,  100.0 };
1626             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1627               Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer] * 0.5
1628                           + 1.0 - cseed[sLayer+iLayer].GetZref(0);
1629               if (max < minmax[1]) {
1630                 minmax[1] = max;
1631               }
1632               Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer] * 0.5
1633                           - 1.0 - cseed[sLayer+iLayer].GetZref(0);
1634               if (min > minmax[0]) {
1635                 minmax[0] = min;
1636               }
1637             }
1638
1639             Bool_t isFake = kFALSE;
1640             if (cl0->GetLabel(0) != cl3->GetLabel(0)) {
1641               isFake = kTRUE;
1642             }
1643             if (cl1->GetLabel(0) != cl3->GetLabel(0)) {
1644               isFake = kTRUE;
1645             }
1646             if (cl2->GetLabel(0) != cl3->GetLabel(0)) {
1647               isFake = kTRUE;
1648             }
1649
1650             if (AliTRDReconstructor::StreamLevel() > 0) {
1651               if ((!isFake) || ((icl3%10) == 0)) {  // Debugging print
1652                 TTreeSRedirector &cstream = *fDebugStreamer;
1653                 cstream << "Seeds0"
1654                         << "isFake=" << isFake
1655                         << "Cl0.="   << cl0
1656                         << "Cl1.="   << cl1
1657                         << "Cl2.="   << cl2
1658                         << "Cl3.="   << cl3
1659                         << "Xref="   << xref
1660                         << "X0="     << xcl[sLayer+0]
1661                         << "X1="     << xcl[sLayer+1]
1662                         << "X2="     << xcl[sLayer+2]
1663                         << "X3="     << xcl[sLayer+3]
1664                         << "Y2exp="  << y2exp
1665                         << "Z2exp="  << z2exp
1666                         << "Chi2R="  << chi2R
1667                         << "Chi2Z="  << chi2Z
1668                         << "Seed0.=" << &cseed[sLayer+0]
1669                         << "Seed1.=" << &cseed[sLayer+1]
1670                         << "Seed2.=" << &cseed[sLayer+2]
1671                         << "Seed3.=" << &cseed[sLayer+3]
1672                         << "Zmin="   << minmax[0]
1673                         << "Zmax="   << minmax[1]
1674                         << "\n";
1675               }
1676             }
1677
1678             ////////////////////////////////////////////////////////////////////////////////////
1679             //
1680             // Fit seeding part
1681             //
1682             ////////////////////////////////////////////////////////////////////////////////////
1683
1684             cl[sLayer+0] = cl0;
1685             cl[sLayer+1] = cl1;
1686             cl[sLayer+2] = cl2;
1687             cl[sLayer+3] = cl3;
1688             Bool_t isOK = kTRUE;
1689
1690             for (Int_t jLayer = 0; jLayer < 4; jLayer++) {
1691
1692               cseed[sLayer+jLayer].SetTilt(hL[sLayer+jLayer]);
1693               cseed[sLayer+jLayer].SetPadLength(padlength[sLayer+jLayer]);
1694               cseed[sLayer+jLayer].SetX0(xcl[sLayer+jLayer]);
1695
1696               for (Int_t iter = 0; iter < 2; iter++) {
1697
1698                 //
1699                 // In iteration 0 we try only one pad-row
1700                 // If quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
1701                 //
1702                 AliTRDseed tseed = cseed[sLayer+jLayer];
1703                 Float_t    roadz = padlength[sLayer+jLayer] * 0.5;
1704                 if (iter > 0) {
1705                   roadz = padlength[sLayer+jLayer];
1706                 }
1707
1708                 Float_t quality = 10000.0;
1709
1710                 for (Int_t iTime = 2; iTime < 20; iTime++) {
1711
1712                   AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1713                   Double_t dxlayer = layer.GetX() - xcl[sLayer+jLayer];
1714                   Double_t zexp    = cl[sLayer+jLayer]->GetZ();
1715
1716                   if (iter > 0) {
1717                     // Try 2 pad-rows in second iteration
1718                     zexp = tseed.GetZref(0) + tseed.GetZref(1) * dxlayer;
1719                     if (zexp > cl[sLayer+jLayer]->GetZ()) {
1720                       zexp = cl[sLayer+jLayer]->GetZ() + padlength[sLayer+jLayer]*0.5;
1721                     }
1722                     if (zexp < cl[sLayer+jLayer]->GetZ()) {
1723                       zexp = cl[sLayer+jLayer]->GetZ() - padlength[sLayer+jLayer]*0.5;
1724                     }
1725                   }
1726
1727                   Double_t yexp  = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
1728                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
1729                   if (index <= 0) {
1730                     continue;
1731                   }
1732                   AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
1733
1734                   tseed.SetIndexes(iTime,index);
1735                   tseed.SetClusters(iTime,cl);  // Register cluster
1736                   tseed.SetX(iTime,dxlayer);    // Register cluster
1737                   tseed.SetY(iTime,cl->GetY()); // Register cluster
1738                   tseed.SetZ(iTime,cl->GetZ()); // Register cluster
1739
1740                 }
1741
1742                 tseed.Update();
1743
1744                 // Count the number of clusters and distortions into quality
1745                 Float_t dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
1746                 Float_t tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle) / 0.1
1747                                  + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0))   / 0.2
1748                                  + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1749                 if ((iter == 0) && tseed.IsOK()) {
1750                   cseed[sLayer+jLayer] = tseed;
1751                   quality              = tquality;
1752                   if (tquality < 5) {
1753                     break;
1754                   }
1755                 }
1756                 if (tseed.IsOK() && (tquality < quality)) {
1757                   cseed[sLayer+jLayer] = tseed;
1758                 }
1759
1760               } // Loop: iter
1761
1762               if (!cseed[sLayer+jLayer].IsOK()) {
1763                 isOK = kFALSE;
1764                 break;
1765               }
1766
1767               cseed[sLayer+jLayer].CookLabels();
1768               cseed[sLayer+jLayer].UpdateUsed();
1769               nusedCl += cseed[sLayer+jLayer].GetNUsed();
1770               if (nusedCl > 25) {
1771                 isOK = kFALSE;
1772                 break;
1773               }
1774
1775             } // Loop: jLayer
1776
1777             if (!isOK) {
1778               continue;
1779             }
1780             nclusters = 0;
1781             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1782               if (cseed[sLayer+iLayer].IsOK()) {
1783                 nclusters += cseed[sLayer+iLayer].GetN2();
1784               }
1785             }
1786
1787             // Iteration 0
1788             rieman.Reset();
1789             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1790               rieman.AddPoint(xcl[sLayer+iLayer]
1791                              ,cseed[sLayer+iLayer].GetYfitR(0)
1792                              ,cseed[sLayer+iLayer].GetZProb()
1793                              ,1
1794                              ,10);
1795             }
1796             rieman.Update();
1797
1798             chi2R = 0.0;
1799             chi2Z = 0.0;
1800
1801             for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1802               cseed[sLayer+iLayer].SetYref(0,rieman.GetYat(xcl[sLayer+iLayer]));
1803               chi2R += (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0))
1804                      * (cseed[sLayer+iLayer].GetYref(0) - cseed[sLayer+iLayer].GetYfitR(0));
1805               cseed[sLayer+iLayer].SetYref(1,rieman.GetDYat(xcl[sLayer+iLayer]));
1806               cseed[sLayer+iLayer].SetZref(0,rieman.GetZat(xcl[sLayer+iLayer]));
1807               chi2Z += (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz())
1808                      * (cseed[sLayer+iLayer].GetZref(0) - cseed[sLayer+iLayer].GetMeanz());
1809               cseed[sLayer+iLayer].SetZref(1,rieman.GetDZat(xcl[sLayer+iLayer]));
1810             }
1811             Double_t curv = rieman.GetC();
1812
1813             //
1814             // Likelihoods
1815             //
1816             Double_t sumda     = TMath::Abs(cseed[sLayer+0].GetYfitR(1) - cseed[sLayer+0].GetYref(1))
1817                                + TMath::Abs(cseed[sLayer+1].GetYfitR(1) - cseed[sLayer+1].GetYref(1))
1818                                + TMath::Abs(cseed[sLayer+2].GetYfitR(1) - cseed[sLayer+2].GetYref(1))
1819                                + TMath::Abs(cseed[sLayer+3].GetYfitR(1) - cseed[sLayer+3].GetYref(1));
1820             Double_t likea     = TMath::Exp(-sumda*10.6);
1821             Double_t likechi2  = 0.0000000001;
1822             if (chi2R < 0.5) {
1823               likechi2 += TMath::Exp(-TMath::Sqrt(chi2R) * 7.73);
1824             }
1825             Double_t likechi2z = TMath::Exp(-chi2Z * 0.088) / TMath::Exp(-chi2Z * 0.019);
1826             Double_t likeN     = TMath::Exp(-(72 - nclusters) * 0.19);
1827             Double_t like      = likea * likechi2 * likechi2z * likeN;
1828             Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetYref(1) - 130.0*curv) * 1.9);
1829             Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].GetZref(1)
1830                                                       - cseed[sLayer+0].GetZref(0) / xcl[sLayer+0]) * 5.9);
1831             Double_t likePrim  = TMath::Max(likePrimY*likePrimZ,0.0005);
1832
1833             seedquality[registered] = like;
1834             seedlayer[registered]   = sLayer;
1835             if (TMath::Log(0.000000000000001 + like) < -15) {
1836               continue;
1837             }
1838             AliTRDseed seedb[6];
1839             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1840               seedb[iLayer] = cseed[iLayer];
1841             }
1842
1843             ////////////////////////////////////////////////////////////////////////////////////
1844             //
1845             // Full track fit part
1846             //
1847             ////////////////////////////////////////////////////////////////////////////////////
1848
1849             Int_t nlayers  = 0;
1850             Int_t nusedf   = 0;
1851             Int_t findable = 0;
1852
1853             //
1854             // Add new layers  - avoid long extrapolation
1855             //
1856             Int_t tLayer[2] = { 0, 0 };
1857             if (sLayer == 2) {
1858               tLayer[0] = 1;
1859               tLayer[1] = 0;
1860             }
1861             if (sLayer == 1) {
1862               tLayer[0] = 5;
1863               tLayer[1] = 0;
1864             }
1865             if (sLayer == 0) {
1866               tLayer[0] = 4;
1867               tLayer[1] = 5;
1868             }
1869
1870             for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1871               Int_t jLayer = tLayer[iLayer]; // Set tracking layer
1872               cseed[jLayer].Reset();
1873               cseed[jLayer].SetTilt(hL[jLayer]);
1874               cseed[jLayer].SetPadLength(padlength[jLayer]);
1875               cseed[jLayer].SetX0(xcl[jLayer]);
1876               // Get pad length and rough cluster
1877               Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].GetYref(0)
1878                                                                       ,cseed[jLayer].GetZref(0)
1879                                                                       ,kRoad2y
1880                                                                       ,kRoad2z);
1881               if (indexdummy <= 0) {
1882                 continue;
1883               }
1884               AliTRDcluster *cldummy = (AliTRDcluster *) GetCluster(indexdummy);
1885               padlength[jLayer]      = TMath::Sqrt(cldummy->GetSigmaZ2() * 12.0);
1886             }
1887             AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1888
1889             for (Int_t iLayer = 0; iLayer < 2; iLayer++) {
1890
1891               Int_t    jLayer = tLayer[iLayer];      // set tracking layer
1892               if ((jLayer == 0) && !(cseed[1].IsOK())) {
1893                 continue;  // break not allowed
1894               }
1895               if ((jLayer == 5) && !(cseed[4].IsOK())) {
1896                 continue;  // break not allowed
1897               }
1898               Float_t  zexp   = cseed[jLayer].GetZref(0);
1899               Double_t zroad  = padlength[jLayer] * 0.5 + 1.0;
1900
1901               for (Int_t iter = 0; iter < 2; iter++) {
1902
1903                 AliTRDseed tseed = cseed[jLayer];
1904                 Float_t quality = 10000.0;
1905
1906                 for (Int_t iTime = 2; iTime < 20; iTime++) {
1907                   AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1908                   Double_t dxlayer = layer.GetX()-xcl[jLayer];
1909                   Double_t yexp    = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer;
1910                   Float_t  yroad   = kRoad1y;
1911                   Int_t    index   = layer.FindNearestCluster(yexp,zexp,yroad,zroad);
1912                   if (index <= 0) {
1913                     continue;
1914                   }
1915                   AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
1916                   tseed.SetIndexes(iTime,index);
1917                   tseed.SetClusters(iTime,cl);  // Register cluster
1918                   tseed.SetX(iTime,dxlayer);    // Register cluster
1919                   tseed.SetY(iTime,cl->GetY()); // Register cluster
1920                   tseed.SetZ(iTime,cl->GetZ()); // Register cluster
1921                 }
1922
1923                 tseed.Update();
1924                 if (tseed.IsOK()) {
1925                   Float_t dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
1926                   Float_t tquality = (18.0 - tseed.GetN2())/2.0 + TMath::Abs(dangle) / 0.1
1927                                    + TMath::Abs(tseed.GetYfit(0) - tseed.GetYref(0)) / 0.2
1928                                    + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0)) / padlength[jLayer];
1929                   if (tquality < quality) {
1930                     cseed[jLayer] = tseed;
1931                     quality       = tquality;
1932                   }
1933                 }
1934
1935                 zroad *= 2.0;
1936
1937               } // Loop: iter
1938
1939               if ( cseed[jLayer].IsOK()) {
1940                 cseed[jLayer].CookLabels();
1941                 cseed[jLayer].UpdateUsed();
1942                 nusedf += cseed[jLayer].GetNUsed();
1943                 AliTRDseed::FitRiemanTilt(cseed,kTRUE);
1944               }
1945
1946             } // Loop: iLayer
1947
1948             // Make copy
1949             AliTRDseed bseed[6];
1950             for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1951               bseed[jLayer] = cseed[jLayer];
1952             }
1953             Float_t lastquality = 10000.0;
1954             Float_t lastchi2    = 10000.0;
1955             Float_t chi2        =  1000.0;
1956
1957             for (Int_t iter = 0; iter < 4; iter++) {
1958
1959               // Sort tracklets according "quality", try to "improve" 4 worst
1960               Float_t sumquality = 0.0;
1961               Float_t squality[6];
1962               Int_t   sortindexes[6];
1963
1964               for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1965                 if (bseed[jLayer].IsOK()) {
1966                   AliTRDseed &tseed = bseed[jLayer];
1967                   Double_t zcor     = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
1968                   Float_t  dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
1969                   Float_t  tquality = (18.0 - tseed.GetN2()) / 2.0 + TMath::Abs(dangle)        / 0.1
1970                                     + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
1971                                     + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0))    / padlength[jLayer];
1972                   squality[jLayer]  = tquality;
1973                 }
1974                 else {
1975                   squality[jLayer] = -1.0;
1976                 }
1977                 sumquality +=squality[jLayer];
1978               }
1979
1980               if ((sumquality >= lastquality) ||
1981                   (chi2       >     lastchi2)) {
1982                 break;
1983               }
1984               lastquality = sumquality;
1985               lastchi2    = chi2;
1986               if (iter > 0) {
1987                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1988                   cseed[jLayer] = bseed[jLayer];
1989                 }
1990               }
1991               TMath::Sort(6,squality,sortindexes,kFALSE);
1992
1993               for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1994
1995                 Int_t bLayer = sortindexes[jLayer];
1996                 AliTRDseed tseed = bseed[bLayer];
1997
1998                 for (Int_t iTime = 2; iTime < 20; iTime++) {
1999
2000                   AliTRDpropagationLayer &layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
2001                   Double_t dxlayer = layer.GetX() - xcl[bLayer];
2002                   Double_t zexp    = tseed.GetZref(0);
2003                   Double_t zcor    = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2004                   Float_t  roadz   = padlength[bLayer] + 1;
2005                   if (TMath::Abs(tseed.GetZProb() - zexp) > 0.5*padlength[bLayer]) {
2006                     roadz = padlength[bLayer] * 0.5;
2007                   }
2008                   if (tseed.GetZfit(1)*tseed.GetZref(1)   < 0.0) {
2009                     roadz = padlength[bLayer] * 0.5;
2010                   }
2011                   if (TMath::Abs(tseed.GetZProb() - zexp) < 0.1*padlength[bLayer]) {
2012                     zexp  = tseed.GetZProb();
2013                     roadz = padlength[bLayer] * 0.5;
2014                   }
2015
2016                   Double_t yexp  = tseed.GetYref(0) + tseed.GetYref(1) * dxlayer - zcor;
2017                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y,roadz);
2018                   if (index <= 0) {
2019                     continue;
2020                   }
2021                   AliTRDcluster *cl = (AliTRDcluster *) GetCluster(index);
2022
2023                   tseed.SetIndexes(iTime,index);
2024                   tseed.SetClusters(iTime,cl);  // Register cluster
2025                   tseed.SetX(iTime,dxlayer);    // Register cluster
2026                   tseed.SetY(iTime,cl->GetY()); // Register cluster
2027                   tseed.SetZ(iTime,cl->GetZ()); // Register cluster
2028
2029                 }
2030
2031                 tseed.Update();
2032                 if (tseed.IsOK()) {
2033                   Float_t  dangle   = tseed.GetYfit(1) - tseed.GetYref(1);
2034                   Double_t zcor     = tseed.GetTilt() * (tseed.GetZProb() - tseed.GetZref(0));
2035                   Float_t  tquality = (18.0 - tseed.GetN2()) / 2.0
2036                                     + TMath::Abs(dangle)     / 0.1
2037                                     + TMath::Abs(tseed.GetYfit(0) - (tseed.GetYref(0) - zcor)) / 0.2
2038                                     + 2.0 * TMath::Abs(tseed.GetMeanz() - tseed.GetZref(0))    / padlength[jLayer];
2039                   if (tquality<squality[bLayer]) {
2040                     bseed[bLayer] = tseed;
2041                   }
2042                 }
2043
2044               } // Loop: jLayer
2045
2046               chi2 = AliTRDseed::FitRiemanTilt(bseed,kTRUE);
2047
2048             } // Loop: iter
2049
2050             nclusters = 0;
2051             nlayers   = 0;
2052             findable  = 0;
2053             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2054               if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) {
2055                 findable++;
2056               }
2057               if (cseed[iLayer].IsOK()) {
2058                 nclusters += cseed[iLayer].GetN2();
2059                 nlayers++;
2060               }
2061             }
2062             if (nlayers < 3) {
2063               continue;
2064             }
2065             rieman.Reset();
2066             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2067               if (cseed[iLayer].IsOK()) {
2068                 rieman.AddPoint(xcl[iLayer]
2069                                ,cseed[iLayer].GetYfitR(0)
2070                                ,cseed[iLayer].GetZProb()
2071                                ,1
2072                                ,10);
2073               }
2074             }
2075             rieman.Update();
2076
2077             chi2RF = 0.0;
2078             chi2ZF = 0.0;
2079             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2080               if (cseed[iLayer].IsOK()) {
2081                 cseed[iLayer].SetYref(0,rieman.GetYat(xcl[iLayer]));
2082                 chi2RF += (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0))
2083                         * (cseed[iLayer].GetYref(0) - cseed[iLayer].GetYfitR(0));
2084                 cseed[iLayer].SetYref(1,rieman.GetDYat(xcl[iLayer]));
2085                 cseed[iLayer].SetZref(0,rieman.GetZat(xcl[iLayer]));
2086                 chi2ZF += (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz())
2087                         * (cseed[iLayer].GetZref(0) - cseed[iLayer].GetMeanz());
2088                 cseed[iLayer].SetZref(1,rieman.GetDZat(xcl[iLayer]));
2089               }
2090             }
2091             chi2RF /= TMath::Max((nlayers - 3.0),1.0);
2092             chi2ZF /= TMath::Max((nlayers - 3.0),1.0);
2093             curv    = rieman.GetC();
2094
2095             Double_t xref2 = (xcl[2] + xcl[3]) * 0.5; // Middle of the chamber
2096             Double_t dzmf  = rieman.GetDZat(xref2);
2097             Double_t zmf   = rieman.GetZat(xref2);
2098
2099             //
2100             // Fit hyperplane
2101             //
2102             Int_t npointsT = 0;
2103             fitterTC.ClearPoints();
2104             fitterT2.ClearPoints();
2105             rieman2.Reset();
2106
2107             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2108
2109               if (!cseed[iLayer].IsOK()) {
2110                 continue;
2111               }
2112
2113               for (Int_t itime = 0; itime < 25; itime++) {
2114
2115                 if (!cseed[iLayer].IsUsable(itime)) {
2116                   continue;
2117                 }
2118                 // X relative to the middle chamber
2119                 Double_t x  = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
2120                 Double_t y  = cseed[iLayer].GetY(itime);
2121                 Double_t z  = cseed[iLayer].GetZ(itime);
2122                 // ExB correction to the correction
2123                 // Tilted rieman
2124                 Double_t uvt[6];
2125                 // Global x
2126                 Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
2127                 Double_t t  = 1.0 / (x2*x2 + y*y);
2128                 uvt[1] = t;                 // t
2129                 uvt[0] = 2.0 * x2 * uvt[1]; // u
2130                 uvt[2] = 2.0 * hL[iLayer] * uvt[1];
2131                 uvt[3] = 2.0 * hL[iLayer] * x * uvt[1];
2132                 uvt[4] = 2.0 * (y + hL[iLayer]*z) * uvt[1];
2133                 Double_t error = 2.0 * 0.2 * uvt[1];
2134                 fitterT2.AddPoint(uvt,uvt[4],error);
2135
2136                 //
2137                 // Constrained rieman
2138                 //
2139                 z = cseed[iLayer].GetZ(itime);
2140                 uvt[0] = 2.0 * x2 * t; // u
2141                 uvt[1] = 2.0 * hL[iLayer] * x2 * uvt[1];
2142                 uvt[2] = 2.0 * (y + hL[iLayer] * (z - GetZ())) * t;
2143                 fitterTC.AddPoint(uvt,uvt[2],error);
2144                 rieman2.AddPoint(x2,y,z,1,10);
2145                 npointsT++;
2146
2147               }
2148
2149             } // Loop: iLayer
2150
2151             rieman2.Update();
2152             fitterTC.Eval();
2153             fitterT2.Eval();
2154             Double_t rpolz0 = fitterT2.GetParameter(3);
2155             Double_t rpolz1 = fitterT2.GetParameter(4);
2156
2157             //
2158             // Linear fitter  - not possible to make boundaries
2159             // Do not accept non possible z and dzdx combinations
2160             //
2161             Bool_t acceptablez = kTRUE;
2162             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2163               if (cseed[iLayer].IsOK()) {
2164                 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2165                 if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > padlength[iLayer] * 0.5 + 1.0) {
2166                   acceptablez = kFALSE;
2167                 }
2168               }
2169             }
2170             if (!acceptablez) {
2171               fitterT2.FixParameter(3,zmf);
2172               fitterT2.FixParameter(4,dzmf);
2173               fitterT2.Eval();
2174               fitterT2.ReleaseParameter(3);
2175               fitterT2.ReleaseParameter(4);
2176               rpolz0 = fitterT2.GetParameter(3);
2177               rpolz1 = fitterT2.GetParameter(4);
2178             }
2179
2180             Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
2181             Double_t chi2TC = fitterTC.GetChisquare() / Float_t(npointsT);
2182             Double_t polz1c = fitterTC.GetParameter(2);
2183             Double_t polz0c = polz1c * xref2;
2184             Double_t aC     =  fitterTC.GetParameter(0);
2185             Double_t bC     =  fitterTC.GetParameter(1);
2186             Double_t cC     =  aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
2187             Double_t aR     =  fitterT2.GetParameter(0);
2188             Double_t bR     =  fitterT2.GetParameter(1);
2189             Double_t dR     =  fitterT2.GetParameter(2);
2190             Double_t cR     =  1.0 + bR*bR - dR*aR;
2191             Double_t dca    =  0.0;
2192             if (cR > 0.0) {
2193               dca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
2194               cR  =  aR / TMath::Sqrt(cR);
2195             }
2196
2197             Double_t chi2ZT2 = 0.0;
2198             Double_t chi2ZTC = 0.0;
2199             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2200               if (cseed[iLayer].IsOK()) {
2201                 Double_t zT2 = rpolz0 + rpolz1 * (xcl[iLayer] - xref2);
2202                 Double_t zTC = polz0c + polz1c * (xcl[iLayer] - xref2);
2203                 chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
2204                 chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
2205               }
2206             }
2207             chi2ZT2 /= TMath::Max((nlayers - 3.0),1.0);
2208             chi2ZTC /= TMath::Max((nlayers - 3.0),1.0);
2209
2210             AliTRDseed::FitRiemanTilt(cseed,kTRUE);
2211             Float_t sumdaf = 0.0;
2212             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2213               if (cseed[iLayer].IsOK()) {
2214                 sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))
2215                                     / cseed[iLayer].GetSigmaY2());
2216               }
2217             }
2218             sumdaf /= Float_t (nlayers - 2.0);
2219
2220             //
2221             // Likelihoods for full track
2222             //
2223             Double_t likezf     = TMath::Exp(-chi2ZF * 0.14);
2224             Double_t likechi2C  = TMath::Exp(-chi2TC * 0.677);
2225             Double_t likechi2TR = TMath::Exp(-chi2TR * 0.78);
2226             Double_t likeaf     = TMath::Exp(-sumdaf * 3.23);
2227             seedquality2[registered] = likezf * likechi2TR * likeaf;
2228
2229             // Still needed ????
2230 //          Bool_t isGold = kFALSE;
2231 //
2232 //          if (nlayers == 6        && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE;   // gold
2233 //          if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE;   // gold
2234 //          if (isGold &&nusedf<10){
2235 //            for (Int_t jLayer=0;jLayer<6;jLayer++){
2236 //              if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
2237 //                seed[index][jLayer].UseClusters();  //sign gold
2238 //            }
2239 //          }
2240
2241             Int_t index0 = 0;
2242             if (!cseed[0].IsOK()) {
2243               index0 = 1;
2244               if (!cseed[1].IsOK()) {
2245                 index0 = 2;
2246               }
2247             }
2248             seedparams[registered][0] = cseed[index0].GetX0();
2249             seedparams[registered][1] = cseed[index0].GetYref(0);
2250             seedparams[registered][2] = cseed[index0].GetZref(0);
2251             seedparams[registered][5] = cR;
2252             seedparams[registered][3] = cseed[index0].GetX0() * cR - TMath::Sin(TMath::ATan(cseed[0].GetYref(1)));
2253             seedparams[registered][4] = cseed[index0].GetZref(1)
2254                                       / TMath::Sqrt(1.0 + cseed[index0].GetYref(1) * cseed[index0].GetYref(1));
2255             seedparams[registered][6] = ns;
2256
2257             Int_t labels[12];
2258             Int_t outlab[24];
2259             Int_t nlab = 0;
2260             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2261               if (!cseed[iLayer].IsOK()) {
2262                 continue;
2263               }
2264               if (cseed[iLayer].GetLabels(0) >= 0) {
2265                 labels[nlab] = cseed[iLayer].GetLabels(0);
2266                 nlab++;
2267               }
2268               if (cseed[iLayer].GetLabels(1) >= 0) {
2269                 labels[nlab] = cseed[iLayer].GetLabels(1);
2270                 nlab++;
2271               }
2272             }
2273             Freq(nlab,labels,outlab,kFALSE);
2274             Int_t label     = outlab[0];
2275             Int_t frequency = outlab[1];
2276             for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2277               cseed[iLayer].SetFreq(frequency);
2278               cseed[iLayer].SetC(cR);
2279               cseed[iLayer].SetCC(cC);
2280               cseed[iLayer].SetChi2(chi2TR);
2281               cseed[iLayer].SetChi2Z(chi2ZF);
2282             }
2283
2284             // Debugging print
2285             if (1 || (!isFake)) {
2286               Float_t zvertex = GetZ();
2287               TTreeSRedirector &cstream = *fDebugStreamer;
2288               if (AliTRDReconstructor::StreamLevel() > 0) {
2289                 cstream << "Seeds1"
2290                         << "isFake="     << isFake
2291                         << "Vertex="     << zvertex
2292                         << "Rieman2.="   << &rieman2
2293                         << "Rieman.="    << &rieman
2294                         << "Xref="       << xref
2295                         << "X0="         << xcl[0]
2296                         << "X1="         << xcl[1]
2297                         << "X2="         << xcl[2]
2298                         << "X3="         << xcl[3]
2299                         << "X4="         << xcl[4]
2300                         << "X5="         << xcl[5]
2301                         << "Chi2R="      << chi2R
2302                         << "Chi2Z="      << chi2Z
2303                         << "Chi2RF="     << chi2RF  // Chi2 of trackletes on full track
2304                         << "Chi2ZF="     << chi2ZF  // Chi2 z on tracklets on full track
2305                         << "Chi2ZT2="    << chi2ZT2 // Chi2 z on tracklets on full track  - rieman tilt
2306                         << "Chi2ZTC="    << chi2ZTC // Chi2 z on tracklets on full track  - rieman tilt const
2307                         << "Chi2TR="     << chi2TR  // Chi2 without vertex constrain
2308                         << "Chi2TC="     << chi2TC  // Chi2 with    vertex constrain
2309                         << "C="          << curv    // Non constrained - no tilt correction
2310                         << "DR="         << dR      // DR parameter    - tilt correction
2311                         << "DCA="        << dca     // DCA             - tilt correction
2312                         << "CR="         << cR      // Non constrained curvature - tilt correction
2313                         << "CC="         << cC      // Constrained curvature
2314                         << "Polz0="      << polz0c
2315                         << "Polz1="      << polz1c
2316                         << "RPolz0="     << rpolz0
2317                         << "RPolz1="     << rpolz1
2318                         << "Ncl="        << nclusters
2319                         << "Nlayers="    << nlayers
2320                         << "NUsedS="     << nusedCl
2321                         << "NUsed="      << nusedf
2322                         << "Findable="   << findable
2323                         << "Like="       << like
2324                         << "LikePrim="   << likePrim
2325                         << "Likechi2C="  << likechi2C
2326                         << "Likechi2TR=" << likechi2TR
2327                         << "Likezf="     << likezf
2328                         << "LikeF="      << seedquality2[registered]
2329                         << "S0.="        << &cseed[0]
2330                         << "S1.="        << &cseed[1]
2331                         << "S2.="        << &cseed[2]
2332                         << "S3.="        << &cseed[3]
2333                         << "S4.="        << &cseed[4]
2334                         << "S5.="        << &cseed[5]
2335                         << "SB0.="       << &seedb[0]
2336                         << "SB1.="       << &seedb[1]
2337                         << "SB2.="       << &seedb[2]
2338                         << "SB3.="       << &seedb[3]
2339                         << "SB4.="       << &seedb[4]
2340                         << "SB5.="       << &seedb[5]
2341                         << "Label="      << label
2342                         << "Freq="       << frequency
2343                         << "sLayer="     << sLayer
2344                         << "\n";
2345               }
2346             }
2347
2348             if (registered<kMaxSeed - 1) {
2349               registered++;
2350               cseed = seed[registered];
2351             }
2352
2353           } // End of loop over layer 1
2354
2355         } // End of loop over layer 0
2356
2357       } // End of loop over layer 3
2358
2359     } // End of loop over seeding time bins
2360
2361     //
2362     // Choose best
2363     //
2364
2365     TMath::Sort(registered,seedquality2,sort,kTRUE);
2366     Bool_t signedseed[kMaxSeed];
2367     for (Int_t i = 0; i < registered; i++) {
2368       signedseed[i] = kFALSE;
2369     }
2370
2371     for (Int_t iter = 0; iter < 5; iter++) {
2372
2373       for (Int_t iseed = 0; iseed < registered; iseed++) {
2374
2375         Int_t index = sort[iseed];
2376         if (signedseed[index]) {
2377           continue;
2378         }
2379         Int_t labelsall[1000];
2380         Int_t nlabelsall = 0;
2381         Int_t naccepted  = 0;;
2382         Int_t sLayer     = seedlayer[index];
2383         Int_t ncl        = 0;
2384         Int_t nused      = 0;
2385         Int_t nlayers    = 0;
2386         Int_t findable   = 0;
2387
2388         for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2389
2390           if (TMath::Abs(seed[index][jLayer].GetYref(0) / xcl[jLayer]) < 0.15) {
2391             findable++;
2392           }
2393           if (seed[index][jLayer].IsOK()) {
2394             seed[index][jLayer].UpdateUsed();
2395             ncl   +=seed[index][jLayer].GetN2();
2396             nused +=seed[index][jLayer].GetNUsed();
2397             nlayers++;
2398             // Cooking label
2399             for (Int_t itime = 0; itime < 25; itime++) {
2400               if (seed[index][jLayer].IsUsable(itime)) {
2401                 naccepted++;
2402                 for (Int_t ilab = 0; ilab < 3; ilab++) {
2403                   Int_t tindex = seed[index][jLayer].GetClusters(itime)->GetLabel(ilab);
2404                   if (tindex >= 0) {
2405                     labelsall[nlabelsall] = tindex;
2406                     nlabelsall++;
2407                   }
2408                 }
2409               }
2410             }
2411           }
2412
2413         }
2414
2415         if (nused > 30) {
2416           continue;
2417         }
2418
2419         if (iter == 0) {
2420           if (nlayers < 6) {
2421             continue;
2422           }
2423           if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2424             continue; // Gold
2425           }
2426         }
2427
2428         if (iter == 1) {
2429           if (nlayers < findable) {
2430             continue;
2431           }
2432           if (TMath::Log(0.000000001+seedquality2[index]) < -4.0) {
2433             continue;
2434           }
2435         }
2436
2437         if (iter == 2) {
2438           if ((nlayers == findable) ||
2439               (nlayers ==        6)) {
2440             continue;
2441           }
2442           if (TMath::Log(0.000000001+seedquality2[index]) < -6.0) {
2443             continue;
2444           }
2445         }
2446
2447         if (iter == 3) {
2448           if (TMath::Log(0.000000001+seedquality2[index]) < -5.0) {
2449             continue;
2450           }
2451         }
2452
2453         if (iter == 4) {
2454           if (TMath::Log(0.000000001+seedquality2[index]) - nused/(nlayers-3.0) < -15.0) {
2455             continue;
2456           }
2457         }
2458
2459         signedseed[index] = kTRUE;
2460
2461         Int_t labels[1000];
2462         Int_t outlab[1000];
2463         Int_t nlab = 0;
2464         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2465           if (seed[index][iLayer].IsOK()) {
2466             if (seed[index][iLayer].GetLabels(0) >= 0) {
2467               labels[nlab] = seed[index][iLayer].GetLabels(0);
2468               nlab++;
2469             }
2470             if (seed[index][iLayer].GetLabels(1) >= 0) {
2471               labels[nlab] = seed[index][iLayer].GetLabels(1);
2472               nlab++;
2473             }
2474           }
2475         }
2476         Freq(nlab,labels,outlab,kFALSE);
2477         Int_t   label     = outlab[0];
2478         Int_t   frequency = outlab[1];
2479         Freq(nlabelsall,labelsall,outlab,kFALSE);
2480         Int_t   label1    = outlab[0];
2481         Int_t   label2    = outlab[2];
2482         Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
2483         Float_t ratio     = Float_t(nused) / Float_t(ncl);
2484         if (ratio < 0.25) {
2485           for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2486             if ((seed[index][jLayer].IsOK()) &&
2487                 (TMath::Abs(seed[index][jLayer].GetYfit(1) - seed[index][jLayer].GetYfit(1)) < 0.2)) {
2488               seed[index][jLayer].UseClusters(); // Sign gold
2489             }
2490           }
2491         }
2492
2493         Int_t eventNr = esd->GetEventNumber();
2494         TTreeSRedirector &cstream = *fDebugStreamer;
2495
2496         //
2497         // Register seed
2498         //
2499         AliTRDtrack *track = RegisterSeed(seed[index],seedparams[index]);
2500         AliTRDtrack  dummy;
2501         if (!track) {
2502           track = &dummy;
2503         }
2504         else {
2505           AliESDtrack esdtrack;
2506           esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
2507           esdtrack.SetLabel(label);
2508           esd->AddTrack(&esdtrack);
2509           TTreeSRedirector &cstream = *fDebugStreamer;
2510           if (AliTRDReconstructor::StreamLevel() > 0) {
2511             cstream << "Tracks"
2512                     << "EventNr="  << eventNr
2513                     << "ESD.="     << &esdtrack
2514                     << "trd.="     << track
2515                     << "trdback.=" << track
2516                     << "\n";
2517           }
2518         }
2519
2520         if (AliTRDReconstructor::StreamLevel() > 0) {
2521           cstream << "Seeds2"
2522                   << "Iter="      << iter
2523                   << "Track.="    << track
2524                   << "Like="      << seedquality[index]
2525                   << "LikeF="     << seedquality2[index]
2526                   << "S0.="       << &seed[index][0]
2527                   << "S1.="       << &seed[index][1]
2528                   << "S2.="       << &seed[index][2]
2529                   << "S3.="       << &seed[index][3]
2530                   << "S4.="       << &seed[index][4]
2531                   << "S5.="       << &seed[index][5]
2532                   << "Label="     << label
2533                   << "Label1="    << label1
2534                   << "Label2="    << label2
2535                   << "FakeRatio=" << fakeratio
2536                   << "Freq="      << frequency
2537                   << "Ncl="       << ncl
2538                   << "Nlayers="   << nlayers
2539                   << "Findable="  << findable
2540                   << "NUsed="     << nused
2541                   << "sLayer="    << sLayer
2542                   << "EventNr="   << eventNr
2543                   << "\n";
2544         }
2545
2546       } // Loop: iseed
2547
2548     } // Loop: iter
2549
2550   } // End of loop over sectors
2551
2552   delete [] pseed;
2553
2554 }
2555
2556 //_____________________________________________________________________________
2557 Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
2558 {
2559   //
2560   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2561   // from the file. The names of the cluster tree and branches
2562   // should match the ones used in AliTRDclusterizer::WriteClusters()
2563   //
2564
2565   Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
2566   TObjArray *clusterArray = new TObjArray(nsize+1000);
2567
2568   TBranch *branch = clusterTree->GetBranch("TRDcluster");
2569   if (!branch) {
2570     AliError("Can't get the branch !");
2571     return 1;
2572   }
2573   branch->SetAddress(&clusterArray);
2574
2575   // Loop through all entries in the tree
2576   Int_t nEntries   = (Int_t) clusterTree->GetEntries();
2577   Int_t nbytes     = 0;
2578   AliTRDcluster *c = 0;
2579   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2580
2581     // Import the tree
2582     nbytes += clusterTree->GetEvent(iEntry);
2583
2584     // Get the number of points in the detector
2585     Int_t nCluster = clusterArray->GetEntriesFast();
2586
2587     // Loop through all TRD digits
2588     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2589       c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
2590       AliTRDcluster *co = c;
2591       array->AddLast(co);
2592       clusterArray->RemoveAt(iCluster);
2593     }
2594
2595   }
2596
2597   delete clusterArray;
2598
2599   return 0;
2600
2601 }
2602
2603 //_____________________________________________________________________________
2604 Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
2605 {
2606   //
2607   // Get track space point with index i
2608   // Origin: C.Cheshkov
2609   //
2610
2611   AliTRDcluster *cl = (AliTRDcluster *) fClusters->UncheckedAt(index);
2612   Int_t  idet     = cl->GetDetector();
2613   Int_t  isector  = fGeom->GetSector(idet);
2614   Int_t  ichamber = fGeom->GetChamber(idet);
2615   Int_t  iplan    = fGeom->GetPlane(idet);
2616   Double_t local[3];
2617   local[0] = GetX(isector,iplan,cl->GetLocalTimeBin());
2618   local[1] = cl->GetY();
2619   local[2] = cl->GetZ();
2620   Double_t global[3];
2621   fGeom->RotateBack(idet,local,global);
2622   p.SetXYZ(global[0],global[1],global[2]);
2623   AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
2624   switch (iplan) {
2625   case 0:
2626     iLayer = AliAlignObj::kTRD1;
2627     break;
2628   case 1:
2629     iLayer = AliAlignObj::kTRD2;
2630     break;
2631   case 2:
2632     iLayer = AliAlignObj::kTRD3;
2633     break;
2634   case 3:
2635     iLayer = AliAlignObj::kTRD4;
2636     break;
2637   case 4:
2638     iLayer = AliAlignObj::kTRD5;
2639     break;
2640   case 5:
2641     iLayer = AliAlignObj::kTRD6;
2642     break;
2643   };
2644   Int_t    modId = isector * fGeom->Ncham() + ichamber;
2645   UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
2646   p.SetVolumeID(volid);
2647
2648   return kTRUE;
2649
2650 }
2651
2652 //_____________________________________________________________________________
2653 void AliTRDtracker::CookLabel(AliKalmanTrack *pt, Float_t wrong) const
2654 {
2655   //
2656   // This cooks a label. Mmmmh, smells good...
2657   //
2658
2659   Int_t label = 123456789;
2660   Int_t index;
2661   Int_t i;
2662   Int_t j;
2663   Int_t ncl   = pt->GetNumberOfClusters();
2664
2665   const Int_t kRange = fTrSec[0]->GetOuterTimeBin() + 1;
2666
2667   Bool_t labelAdded;
2668
2669   Int_t **s = new Int_t* [kRange];
2670   for (i = 0; i < kRange; i++) {
2671     s[i] = new Int_t[2];
2672   }
2673   for (i = 0; i < kRange; i++) {
2674     s[i][0] = -1;
2675     s[i][1] =  0;
2676   }
2677
2678   Int_t t0;
2679   Int_t t1;
2680   Int_t t2;
2681
2682   for (i = 0; i < ncl; i++) {
2683     index = pt->GetClusterIndex(i);
2684     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2685     t0=c->GetLabel(0);
2686     t1=c->GetLabel(1);
2687     t2=c->GetLabel(2);
2688   }
2689
2690   for (i = 0; i < ncl; i++) {
2691     index = pt->GetClusterIndex(i);
2692     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2693     for (Int_t k = 0; k < 3; k++) {
2694       label      = c->GetLabel(k);
2695       labelAdded = kFALSE;
2696       j = 0;
2697       if (label >= 0) {
2698         while ((!labelAdded) && (j < kRange)) {
2699           if ((s[j][0] == label) ||
2700               (s[j][1] ==     0)) {
2701             s[j][0] = label;
2702             s[j][1] = s[j][1] + 1;
2703             labelAdded = kTRUE;
2704           }
2705           j++;
2706         }
2707       }
2708     }
2709   }
2710
2711   Int_t max = 0;
2712   label = -123456789;
2713
2714   for (i = 0; i < kRange; i++) {
2715     if (s[i][1] > max) {
2716       max   = s[i][1];
2717       label = s[i][0];
2718     }
2719   }
2720
2721   for (i = 0; i < kRange; i++) {
2722     delete []s[i];
2723   }
2724
2725   delete []s;
2726
2727   if ((1.0 - Float_t(max)/ncl) > wrong) {
2728     label = -label;
2729   }
2730
2731   pt->SetLabel(label);
2732
2733 }
2734
2735 //_____________________________________________________________________________
2736 void AliTRDtracker::UseClusters(const AliKalmanTrack *t, Int_t from) const
2737 {
2738   //
2739   // Use clusters, but don't abuse them!
2740   //
2741
2742   const Float_t kmaxchi2 = 18;
2743   const Float_t kmincl   = 10;
2744
2745   AliTRDtrack *track = (AliTRDtrack *) t;
2746
2747   Int_t ncl = t->GetNumberOfClusters();
2748   for (Int_t i = from; i < ncl; i++) {
2749     Int_t index  = t->GetClusterIndex(i);
2750     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(index);
2751     Int_t iplane = fGeom->GetPlane(c->GetDetector());
2752     if (track->fTracklets[iplane].GetChi2() > kmaxchi2) {
2753       continue;
2754     }
2755     if (track->fTracklets[iplane].GetN()    <   kmincl) {
2756       continue;
2757     }
2758     if (!(c->IsUsed())) {
2759       c->Use();
2760     }
2761   }
2762
2763 }
2764
2765 //_____________________________________________________________________________
2766 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
2767 {
2768   //
2769   // Parametrised "expected" error of the cluster reconstruction in Y
2770   //
2771
2772   Double_t s = 0.08 * 0.08;
2773   return s;
2774
2775 }
2776
2777 //_____________________________________________________________________________
2778 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
2779 {
2780   //
2781   // Parametrised "expected" error of the cluster reconstruction in Z
2782   //
2783
2784   Double_t s = 9.0 * 9.0 / 12.0;
2785   return s;
2786
2787 }
2788
2789 //_____________________________________________________________________________
2790 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
2791 {
2792   //
2793   // Returns radial position which corresponds to time bin <localTB>
2794   // in tracking sector <sector> and plane <plane>
2795   //
2796
2797   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
2798   Int_t pl    = fTrSec[sector]->GetLayerNumber(index);
2799
2800   return fTrSec[sector]->GetLayer(pl)->GetX();
2801
2802 }
2803
2804 //_____________________________________________________________________________
2805 AliTRDtracker::AliTRDpropagationLayer
2806              ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
2807                                     , Double_t radLength, Int_t tbIndex, Int_t plane)
2808   :fN(0)
2809   ,fSec(0)
2810   ,fClusters(NULL)
2811   ,fIndex(NULL)
2812   ,fX(x)
2813   ,fdX(dx)
2814   ,fRho(rho)
2815   ,fX0(radLength)
2816   ,fTimeBinIndex(tbIndex)
2817   ,fPlane(plane)
2818   ,fYmax(0)
2819   ,fYmaxSensitive(0)
2820   ,fHole(kFALSE)
2821   ,fHoleZc(0)
2822   ,fHoleZmax(0)
2823   ,fHoleYc(0)
2824   ,fHoleYmax(0)
2825   ,fHoleRho(0)
2826   ,fHoleX0(0)
2827 {
2828   //
2829   // AliTRDpropagationLayer constructor
2830   //
2831
2832   for (Int_t i = 0; i < (Int_t) kZones; i++) {
2833     fZc[i]   = 0;
2834     fZmax[i] = 0;
2835   }
2836
2837   if (fTimeBinIndex >= 0) {
2838     fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
2839     fIndex    = new UInt_t[kMaxClusterPerTimeBin];
2840   }
2841
2842   for (Int_t i = 0; i < 5; i++) {
2843     fIsHole[i] = kFALSE;
2844   }
2845
2846 }
2847
2848 //_____________________________________________________________________________
2849 void AliTRDtracker::AliTRDpropagationLayer
2850                   ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
2851                           , Double_t radLength, Double_t Yc, Double_t Zc)
2852 {
2853   //
2854   // Sets hole in the layer
2855   //
2856
2857   fHole     = kTRUE;
2858   fHoleZc   = Zc;
2859   fHoleZmax = Zmax;
2860   fHoleYc   = Yc;
2861   fHoleYmax = Ymax;
2862   fHoleRho  = rho;
2863   fHoleX0   = radLength;
2864
2865 }
2866
2867 //_____________________________________________________________________________
2868 AliTRDtracker::AliTRDtrackingSector
2869              ::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs)
2870   :fN(0)
2871   ,fGeom(geo)
2872   ,fGeomSector(gs)
2873 {
2874   //
2875   // AliTRDtrackingSector Constructor
2876   //
2877
2878   AliTRDpadPlane         *padPlane = 0;
2879   AliTRDpropagationLayer *ppl      = 0;
2880
2881   // Get holes description from geometry
2882   Bool_t holes[AliTRDgeometry::kNcham];
2883   for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
2884     holes[icham] = fGeom->IsHole(0,icham,gs);
2885   }
2886
2887   for (UInt_t i = 0; i < kMaxTimeBinIndex; i++) {
2888     fTimeBinIndex[i] = -1;
2889   }
2890
2891   Double_t x;
2892   Double_t dx;
2893   Double_t rho;
2894   Double_t radLength;
2895
2896   // Add layers for each of the planes
2897   Double_t dxAmp   = (Double_t) fGeom->CamHght(); // Amplification region
2898   //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2899
2900   const Int_t kNchambers  = AliTRDgeometry::Ncham();
2901   Int_t    tbIndex;
2902   Double_t ymax           = 0;
2903   Double_t ymaxsensitive  = 0;
2904   Double_t *zc            = new Double_t[kNchambers];
2905   Double_t *zmax          = new Double_t[kNchambers];
2906   Double_t *zmaxsensitive = new Double_t[kNchambers];
2907
2908   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
2909   if (!commonParam) {
2910     AliErrorGeneral("AliTRDtrackingSector::Ctor"
2911                    ,"Could not get common parameters\n");
2912     return;
2913   }
2914
2915   for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2916
2917     ymax          = fGeom->GetChamberWidth(plane) / 2.0;
2918     padPlane      = commonParam->GetPadPlane(plane,0);
2919     ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
2920
2921     for (Int_t ch = 0; ch < kNchambers; ch++) {
2922       zmax[ch]          = fGeom->GetChamberLength(plane,ch) / 2.0;
2923       Float_t pad       = padPlane->GetRowSize(1);
2924       Float_t row0      = commonParam->GetRow0(plane,ch,0);
2925       Int_t   nPads     = commonParam->GetRowMax(plane,ch,0);
2926       zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0;
2927       zc[ch]            = -(pad * nPads) / 2.0 + row0;
2928     }
2929
2930     dx        = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
2931               / AliTRDcalibDB::Instance()->GetSamplingFrequency();
2932     rho       = 0.00295 * 0.85; //????
2933     radLength = 11.0;
2934
2935     Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
2936     //Double_t xbottom = x0 - dxDrift;
2937     //Double_t xtop    = x0 + dxAmp;
2938
2939     Int_t nTimeBins =  AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
2940     for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
2941
2942       Double_t xlayer = iTime * dx - dxAmp;
2943       //if (xlayer<0) xlayer = dxAmp / 2.0;
2944       x = x0 - xlayer;
2945
2946       tbIndex = CookTimeBinIndex(plane,iTime);
2947       ppl     = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,plane);
2948       ppl->SetYmax(ymax,ymaxsensitive);
2949       ppl->SetZ(zc,zmax,zmaxsensitive);
2950       ppl->SetHoles(holes);
2951       InsertLayer(ppl);
2952
2953     }
2954
2955   }
2956
2957   MapTimeBinLayers();
2958
2959   delete [] zc;
2960   delete [] zmax;
2961   delete [] zmaxsensitive;
2962
2963 }
2964
2965 //_____________________________________________________________________________
2966 AliTRDtracker::AliTRDtrackingSector
2967              ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
2968   :fN(0)
2969   ,fGeom(0)
2970   ,fGeomSector(0)
2971 {
2972   //
2973   // Copy constructor
2974   //
2975
2976 }
2977
2978 //_____________________________________________________________________________
2979 Int_t  AliTRDtracker::AliTRDtrackingSector
2980                     ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
2981 {
2982   //
2983   // depending on the digitization parameters calculates "global"
2984   // time bin index for timebin <localTB> in plane <plane>
2985   //
2986   //
2987
2988   Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
2989   Int_t gtb        = (plane+1) * tbPerPlane - localTB - 1;
2990   if (localTB < 0) {
2991     return -1;
2992   }
2993   if (gtb     < 0) {
2994     return -1;
2995   }
2996
2997   return gtb;
2998
2999 }
3000
3001 //_____________________________________________________________________________
3002 void AliTRDtracker::AliTRDtrackingSector
3003                   ::MapTimeBinLayers()
3004 {
3005   //
3006   // For all sensitive time bins sets corresponding layer index
3007   // in the array fTimeBins
3008   //
3009
3010   Int_t index;
3011
3012   for (Int_t i = 0; i < fN; i++) {
3013
3014     index = fLayers[i]->GetTimeBinIndex();
3015
3016     if (index < 0) {
3017       continue;
3018     }
3019     if (index >= (Int_t) kMaxTimeBinIndex) {
3020       //AliWarning(Form("Index %d exceeds allowed maximum of %d!\n"
3021       //               ,index,kMaxTimeBinIndex-1));
3022       continue;
3023     }
3024
3025     fTimeBinIndex[index] = i;
3026
3027   }
3028
3029 }
3030
3031 //_____________________________________________________________________________
3032 Int_t AliTRDtracker::AliTRDtrackingSector
3033                    ::GetLayerNumber(Double_t x) const
3034 {
3035   //
3036   // Returns the number of time bin which in radial position is closest to <x>
3037   //
3038
3039   if (x >= fLayers[fN-1]->GetX()) {
3040     return fN - 1;
3041   }
3042   if (x <= fLayers[   0]->GetX()) {
3043     return 0;
3044   }
3045
3046   Int_t b = 0;
3047   Int_t e = fN - 1;
3048   Int_t m = (b + e) / 2;
3049
3050   for ( ; b < e; m = (b + e) / 2) {
3051     if (x > fLayers[m]->GetX()) {
3052       b = m + 1;
3053     }
3054     else {
3055       e = m;
3056     }
3057   }
3058
3059   if (TMath::Abs(x - fLayers[m]->GetX()) > TMath::Abs(x - fLayers[m+1]->GetX())) {
3060     return m + 1;
3061   }
3062   else {
3063     return m;
3064   }
3065
3066 }
3067
3068 //_____________________________________________________________________________
3069 Int_t AliTRDtracker::AliTRDtrackingSector
3070                    ::GetInnerTimeBin() const
3071 {
3072   //
3073   // Returns number of the innermost SENSITIVE propagation layer
3074   //
3075
3076   return GetLayerNumber(0);
3077
3078 }
3079
3080 //_____________________________________________________________________________
3081 Int_t AliTRDtracker::AliTRDtrackingSector
3082                    ::GetOuterTimeBin() const
3083 {
3084   //
3085   // Returns number of the outermost SENSITIVE time bin
3086   //
3087
3088   return GetLayerNumber(GetNumberOfTimeBins() - 1);
3089
3090 }
3091
3092 //_____________________________________________________________________________
3093 Int_t AliTRDtracker::AliTRDtrackingSector
3094                    ::GetNumberOfTimeBins() const
3095 {
3096   //
3097   // Returns number of SENSITIVE time bins
3098   //
3099
3100   Int_t tb;
3101   Int_t layer;
3102
3103   for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
3104     layer = GetLayerNumber(tb);
3105     if (layer >= 0) {
3106       break;
3107     }
3108   }
3109
3110   return tb + 1;
3111
3112 }
3113
3114 //_____________________________________________________________________________
3115 void AliTRDtracker::AliTRDtrackingSector
3116                   ::InsertLayer(AliTRDpropagationLayer *pl)
3117 {
3118   //
3119   // Insert layer <pl> in fLayers array.
3120   // Layers are sorted according to X coordinate.
3121   //
3122
3123   if (fN == ((Int_t) kMaxLayersPerSector)) {
3124     //AliWarning("Too many layers !\n");
3125     return;
3126   }
3127
3128   if (fN == 0) {
3129     fLayers[fN++] = pl;
3130     return;
3131   }
3132
3133   Int_t i = Find(pl->GetX());
3134
3135   memmove(fLayers+i+1,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3136
3137   fLayers[i] = pl;
3138   fN++;
3139
3140 }
3141
3142 //_____________________________________________________________________________
3143 Int_t AliTRDtracker::AliTRDtrackingSector
3144                    ::Find(Double_t x) const
3145 {
3146   //
3147   // Returns index of the propagation layer nearest to X
3148   //
3149
3150   if (x <= fLayers[0]->GetX()) {
3151     return 0;
3152   }
3153
3154   if (x > fLayers[fN-1]->GetX()) {
3155     return fN;
3156   }
3157
3158   Int_t b = 0;
3159   Int_t e = fN-1;
3160   Int_t m = (b + e) / 2;
3161
3162   for (; b < e; m = (b + e) / 2) {
3163     if (x > fLayers[m]->GetX()) {
3164       b = m + 1;
3165     }
3166     else {
3167       e = m;
3168     }
3169   }
3170
3171   return m;
3172
3173 }
3174
3175 //_____________________________________________________________________________
3176 void AliTRDtracker::AliTRDpropagationLayer
3177                   ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
3178 {
3179   //
3180   // set centers and the width of sectors
3181   //
3182
3183   for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3184     fZc[icham]            = center[icham];
3185     fZmax[icham]          = w[icham];
3186     fZmaxSensitive[icham] = wsensitive[icham];
3187   }
3188
3189 }
3190
3191 //_____________________________________________________________________________
3192 void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3193 {
3194   //
3195   // set centers and the width of sectors
3196   //
3197
3198   fHole = kFALSE;
3199
3200   for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3201     fIsHole[icham] = holes[icham];
3202     if (holes[icham]) {
3203       fHole = kTRUE;
3204     }
3205   }
3206
3207 }
3208
3209 //_____________________________________________________________________________
3210 void AliTRDtracker::AliTRDpropagationLayer
3211                   ::InsertCluster(AliTRDcluster *c, UInt_t index)
3212 {
3213   //
3214   // Insert cluster in cluster array.
3215   // Clusters are sorted according to Y coordinate.
3216   //
3217
3218   if (fTimeBinIndex < 0) {
3219     //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
3220     return;
3221   }
3222
3223   if (fN == (Int_t) kMaxClusterPerTimeBin) {
3224     //AliWarning("Too many clusters !\n");
3225     return;
3226   }
3227
3228   if (fN == 0) {
3229     fIndex[0]       = index;
3230     fClusters[fN++] = c;
3231     return;
3232   }
3233
3234   Int_t i = Find(c->GetY());
3235   memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3236   memmove(fIndex   +i+1,fIndex   +i,(fN-i)*sizeof(UInt_t));
3237   fIndex[i]    = index;
3238   fClusters[i] = c;
3239   fN++;
3240
3241 }
3242
3243 //_____________________________________________________________________________
3244 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
3245 {
3246   //
3247   // Returns index of the cluster nearest in Y
3248   //
3249
3250   if (fN <= 0) {
3251     return 0;
3252   }
3253   if (y <= fClusters[0]->GetY()) {
3254     return 0;
3255   }
3256   if (y >  fClusters[fN-1]->GetY()) {
3257     return fN;
3258   }
3259
3260   Int_t b = 0;
3261   Int_t e = fN - 1;
3262   Int_t m = (b + e) / 2;
3263
3264   for ( ; b < e; m = (b + e) / 2) {
3265     if (y > fClusters[m]->GetY()) {
3266       b = m + 1;
3267     }
3268     else {
3269       e = m;
3270     }
3271   }
3272
3273   return m;
3274
3275 }
3276
3277 //_____________________________________________________________________________
3278 Int_t AliTRDtracker::AliTRDpropagationLayer
3279                    ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
3280                                       , Float_t maxroadz) const
3281 {
3282   //
3283   // Returns index of the cluster nearest to the given y,z
3284   //
3285
3286   Int_t   index   = -1;
3287   Int_t   maxn    = fN;
3288   Float_t mindist = maxroad;
3289
3290   for (Int_t i = Find(y-maxroad); i < maxn; i++) {
3291     AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
3292     Float_t ycl = c->GetY();
3293     if (ycl > (y + maxroad)) {
3294       break;
3295     }
3296     if (TMath::Abs(c->GetZ() - z) > maxroadz) {
3297       continue;
3298     }
3299     if (TMath::Abs(ycl - y)       < mindist) {
3300       mindist = TMath::Abs(ycl - y);
3301       index   = fIndex[i];
3302     }
3303   }
3304
3305   return index;
3306
3307 }
3308
3309 //_____________________________________________________________________________
3310 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
3311 {
3312   //
3313   // Returns correction factor for tilted pads geometry
3314   //
3315
3316   Int_t    det   = c->GetDetector();
3317   Int_t    plane = fGeom->GetPlane(det);
3318   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
3319   Double_t h01   = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
3320
3321   if (fNoTilt) {
3322     h01 = 0;
3323   }
3324
3325   return h01;
3326
3327 }
3328
3329 //_____________________________________________________________________________
3330 void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
3331 {
3332   //
3333   // This is setting fdEdxPlane and fTimBinPlane
3334   // Sums up the charge in each plane for track TRDtrack and also get the
3335   // Time bin for Max. Cluster
3336   // Prashant Shukla (shukla@physi.uni-heidelberg.de)
3337   //
3338
3339   Double_t  clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
3340   Double_t  maxclscharge[AliESDtrack::kNPlane];
3341   Int_t     nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
3342   Int_t     timebin[AliESDtrack::kNPlane];
3343
3344   // Initialization of cluster charge per plane.
3345   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
3346     for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
3347       clscharge[iPlane][iSlice] = 0.0;
3348       nCluster[iPlane][iSlice]  = 0;
3349     }
3350   }
3351
3352   // Initialization of cluster charge per plane.
3353   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
3354     timebin[iPlane]      =  -1;
3355     maxclscharge[iPlane] = 0.0;
3356   }
3357
3358   // Loop through all clusters associated to track TRDtrack
3359   Int_t nClus = TRDtrack.GetNumberOfClusters();  // from Kalmantrack
3360   for (Int_t iClus = 0; iClus < nClus; iClus++) {
3361     Double_t charge = TRDtrack.GetClusterdQdl(iClus);
3362     Int_t    index  = TRDtrack.GetClusterIndex(iClus);
3363     AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index);
3364     if (!pTRDcluster) {
3365       continue;
3366     }
3367     Int_t    tb     = pTRDcluster->GetLocalTimeBin();
3368     if (!tb) {
3369       continue;
3370     }
3371     Int_t detector  = pTRDcluster->GetDetector();
3372     Int_t iPlane    = fGeom->GetPlane(detector);
3373     Int_t iSlice    = tb * AliESDtrack::kNSlice / AliTRDtrack::kNtimeBins;
3374     clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge;
3375     if (charge > maxclscharge[iPlane]) {
3376       maxclscharge[iPlane] = charge;
3377       timebin[iPlane]      = tb;
3378     }
3379     nCluster[iPlane][iSlice]++;
3380   } // End of loop over cluster
3381
3382   // Setting the fdEdxPlane and fTimBinPlane variabales
3383   Double_t totalCharge = 0.0;
3384
3385   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
3386     for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
3387       if (nCluster[iPlane][iSlice]) {
3388         clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
3389       }
3390       TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice);
3391       totalCharge = totalCharge+clscharge[iPlane][iSlice];
3392     }
3393     TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);
3394   }
3395
3396   // Still needed ????
3397   //  Int_t i;
3398   //  Int_t nc=TRDtrack.GetNumberOfClusters();
3399   //  Float_t dedx=0;
3400   //  for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
3401   //  dedx /= nc;
3402   //  for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3403   //    TRDtrack.SetPIDsignals(dedx, iPlane);
3404   //    TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
3405   //  }
3406
3407 }
3408
3409 //_____________________________________________________________________________
3410 Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
3411                                 , AliTRDtrack *track
3412                                 , Int_t *clusters, AliTRDtracklet &tracklet)
3413 {
3414   //
3415   //
3416   // Try to find nearest clusters to the track in timebins from t0 to t1
3417   //
3418   // Correction coeficients - depend on TRD parameters - to be changed accordingly
3419   //
3420
3421   Double_t x[100];
3422   Double_t yt[100];
3423   Double_t zt[100];
3424   Double_t xmean = 0.0;       // Reference x
3425   Double_t dz[10][100];
3426   Double_t dy[10][100];
3427   Float_t  zmean[100];
3428   Float_t  nmean[100];
3429   Int_t    clfound = 0;
3430   Int_t    indexes[10][100];  // Indexes of the clusters in the road
3431   Int_t    best[10][100];     // Index of best matching cluster
3432   AliTRDcluster *cl[10][100]; // Pointers to the clusters in the road
3433
3434   for (Int_t it = 0; it < 100; it++) {
3435     x[it]        = 0.0;
3436     yt[it]       = 0.0;
3437     zt[it]       = 0.0;
3438     clusters[it] = -2;
3439     zmean[it]    = 0.0;
3440     nmean[it]    = 0.0;
3441     for (Int_t ih = 0; ih < 10;ih++) {
3442       indexes[ih][it] = -2;   // Reset indexes1
3443       cl[ih][it]      =  0;
3444       dz[ih][it]      = -100.0;
3445       dy[ih][it]      = -100.0;
3446       best[ih][it]    =  0;
3447     }
3448   }
3449
3450   Double_t x0        = track->GetX();
3451   Double_t sigmaz    = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3452   Int_t    nall      = 0;
3453   Int_t    nfound    = 0;
3454   Double_t h01       = 0.0;
3455   Int_t    plane     = -1;
3456   Int_t    detector  = -1;
3457   Float_t  padlength = 0.0;
3458   AliTRDtrack track2(* track);
3459   Float_t  snpy      = track->GetSnp();
3460   Float_t  tany      = TMath::Sqrt(snpy*snpy / (1.0 - snpy*snpy));
3461   if (snpy < 0.0) {
3462     tany *= -1.0;
3463   }
3464
3465   Double_t sy2       = ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
3466   Double_t sz2       = ExpectedSigmaZ2(x0,track->GetTgl());
3467   Double_t road      = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3468   if (road > 6.0) {
3469     road = 6.0;
3470   }
3471
3472   for (Int_t it = 0; it < t1-t0; it++) {
3473
3474     Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };
3475     AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(it+t0));
3476     if (timeBin == 0) {
3477       continue; // No indexes1
3478     }
3479
3480     Int_t    maxn = timeBin;
3481     x[it]  = timeBin.GetX();
3482     track2.PropagateTo(x[it]);
3483     yt[it] = track2.GetY();
3484     zt[it] = track2.GetZ();
3485
3486     Double_t y    = yt[it];
3487     Double_t z    = zt[it];
3488     Double_t chi2 = 1000000.0;
3489     nall++;
3490
3491     //
3492     // Find 2 nearest cluster at given time bin
3493     //
3494     for (Int_t i = timeBin.Find(y - road); i < maxn; i++) {
3495
3496       AliTRDcluster *c = (AliTRDcluster *) (timeBin[i]);
3497       h01 = GetTiltFactor(c);
3498       if (plane < 0) {
3499         Int_t det = c->GetDetector();
3500         plane     = fGeom->GetPlane(det);
3501         padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3502       }
3503       //if (c->GetLocalTimeBin()==0) continue;
3504       if (c->GetY() > (y + road)) {
3505         break;
3506       }
3507       if ((c->GetZ() - z) * (c->GetZ() - z) > (12.0 * sz2)) {
3508         continue;
3509       }
3510
3511       Double_t dist = TMath::Abs(c->GetZ() - z);
3512       if (dist > (0.5 * padlength + 6.0 * sigmaz)) {
3513         continue;   // 6 sigma boundary cut
3514       }
3515       Double_t cost = 0.0;
3516       // Sigma boundary cost function
3517       if (dist> (0.5 * padlength - sigmaz)){
3518         cost =  (dist - 0.5*padlength) / (2.0 * sigmaz);
3519         if (cost > -1) {
3520           cost = (cost + 1.0) * (cost + 1.0);
3521         }
3522         else {
3523           cost = 0.0;
3524         }
3525       }
3526       //Int_t label = TMath::Abs(track->GetLabel());
3527       //if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3528       chi2 = track2.GetPredictedChi2(c,h01) + cost;
3529       clfound++;
3530
3531       if (chi2 > maxChi2[1]) {
3532         continue;
3533       }
3534       detector = c->GetDetector();
3535       // Store the clusters in the road
3536       for (Int_t ih = 2; ih < 9; ih++) {
3537         if (cl[ih][it] == 0) {
3538           cl[ih][it]      = c;
3539           indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
3540           break;
3541         }
3542       }
3543
3544       if (chi2 < maxChi2[0]) {
3545         maxChi2[1]     = maxChi2[0];
3546         maxChi2[0]     = chi2;
3547         indexes[1][it] = indexes[0][it];
3548         cl[1][it]      = cl[0][it];
3549         indexes[0][it] = timeBin.GetIndex(i);
3550         cl[0][it]      = c;
3551         continue;
3552       }
3553       maxChi2[1]     = chi2;
3554       cl[1][it]      = c;
3555       indexes[1][it] = timeBin.GetIndex(i);
3556
3557     }
3558
3559     if (cl[0][it]) {
3560       nfound++;
3561       xmean += x[it];
3562     }
3563
3564   }
3565
3566   if (nfound < 4) {
3567     return 0;
3568   }
3569   xmean /= Float_t(nfound);  // Middle x
3570   track2.PropagateTo(xmean); // Propagate track to the center
3571
3572   //
3573   // Choose one of the variants
3574   //
3575   Int_t    changes[10];
3576   Float_t  sumz   = 0.0;
3577   Float_t  sum    = 0.0;
3578   Double_t sumdy  = 0.0;
3579   Double_t sumdy2 = 0.0;
3580   Double_t sumx   = 0.0;
3581   Double_t sumxy  = 0.0;
3582   Double_t sumx2  = 0.0;
3583   Double_t mpads  = 0.0;
3584
3585   Int_t    ngood[10];
3586   Int_t    nbad[10];
3587
3588   Double_t meanz[10];
3589   Double_t moffset[10];    // Mean offset
3590   Double_t mean[10];       // Mean value
3591   Double_t angle[10];      // Angle
3592
3593   Double_t smoffset[10];   // Sigma of mean offset
3594   Double_t smean[10];      // Sigma of mean value
3595   Double_t sangle[10];     // Sigma of angle
3596   Double_t smeanangle[10]; // Correlation
3597
3598   Double_t sigmas[10];
3599   Double_t tchi2s[10];     // Chi2s for tracklet
3600
3601   for (Int_t it = 0; it < 10; it++) {
3602
3603     ngood[it]      = 0;
3604     nbad[it]       = 0;
3605
3606     meanz[it]      = 0.0;
3607     moffset[it]    = 0.0;    // Mean offset
3608     mean[it]       = 0.0;    // Mean value
3609     angle[it]      = 0.0;    // Angle
3610
3611     smoffset[it]   = 1.0e10; // Sigma of mean offset
3612     smean[it]      = 1.0e10; // Sigma of mean value
3613     sangle[it]     = 1.0e10; // Sigma of angle
3614     smeanangle[it] = 0.0;    // Correlation
3615
3616     sigmas[it]     = 1.0e10;
3617     tchi2s[it]     = 1.0e10; // Chi2s for tracklet
3618
3619   }
3620
3621   //
3622   // Calculate zmean
3623   //
3624   for (Int_t it = 0; it < t1 - t0; it++) {
3625     if (!cl[0][it]) {
3626       continue;
3627     }
3628     for (Int_t dt = -3; dt <= 3; dt++) {
3629       if (it+dt <     0) {
3630         continue;
3631       }
3632       if (it+dt > t1-t0) {
3633         continue;
3634       }
3635       if (!cl[0][it+dt]) {
3636         continue;
3637       }
3638       zmean[it] += cl[0][it+dt]->GetZ();
3639       nmean[it] += 1.0;
3640     }
3641     zmean[it] /= nmean[it];
3642   }
3643
3644   for (Int_t it = 0; it < t1 - t0; it++) {
3645
3646     best[0][it] = 0;
3647
3648     for (Int_t ih = 0; ih < 10; ih++) {
3649       dz[ih][it] = -100.0;
3650       dy[ih][it] = -100.0;
3651       if (!cl[ih][it]) {
3652         continue;
3653       }
3654       Double_t xcluster = cl[ih][it]->GetX();
3655       Double_t ytrack;
3656       Double_t ztrack;
3657       track2.GetProlongation(xcluster,ytrack,ztrack );
3658       dz[ih][it] = cl[ih][it]->GetZ()- ztrack;                   // Calculate distance from track in z
3659       dy[ih][it] = cl[ih][it]->GetY() + dz[ih][it]*h01 - ytrack; // and in y
3660     }
3661
3662     // Minimize changes
3663     if (!cl[0][it]) {
3664       continue;
3665     }
3666     if ((TMath::Abs(cl[0][it]->GetZ()-zmean[it]) > padlength * 0.8) &&
3667         (cl[1][it])) {
3668       if (TMath::Abs(cl[1][it]->GetZ()-zmean[it]) < padlength * 0.5) {
3669         best[0][it] = 1;
3670       }
3671     }
3672
3673   }
3674
3675   //
3676   // Iterative choice of "best path"
3677   //
3678   Int_t label    = TMath::Abs(track->GetLabel());
3679   Int_t bestiter = 0;
3680
3681   for (Int_t iter = 0; iter < 9; iter++) {
3682
3683     changes[iter] = 0;
3684     sumz          = 0;
3685     sum           = 0;
3686     sumdy         = 0;
3687     sumdy2        = 0;
3688     sumx          = 0;
3689     sumx2         = 0;
3690     sumxy         = 0;
3691     mpads         = 0;
3692     ngood[iter]   = 0;
3693     nbad[iter]    = 0;
3694
3695     // Linear fit
3696     for (Int_t it = 0; it < t1 - t0; it++) {
3697
3698       if (!cl[best[iter][it]][it]) {
3699         continue;
3700       }
3701
3702       // Calculates pad-row changes
3703       Double_t zbefore = cl[best[iter][it]][it]->GetZ();
3704       Double_t zafter  = cl[best[iter][it]][it]->GetZ();
3705       for (Int_t itd = it - 1; itd >= 0; itd--) {
3706         if (cl[best[iter][itd]][itd]) {
3707           zbefore = cl[best[iter][itd]][itd]->GetZ();
3708           break;
3709         }
3710       }
3711       for (Int_t itd = it + 1; itd < t1 - t0; itd++) {
3712         if (cl[best[iter][itd]][itd]) {
3713           zafter  = cl[best[iter][itd]][itd]->GetZ();
3714           break;
3715         }
3716       }
3717       if ((TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore) > 0.1) &&
3718           (TMath::Abs(cl[best[iter][it]][it]->GetZ()- zafter) > 0.1)) {
3719         changes[iter]++;
3720       }
3721
3722       Double_t dx = x[it]-xmean; // Distance to reference x
3723       sumz   += cl[best[iter][it]][it]->GetZ();
3724       sum++;
3725       sumdy  += dy[best[iter][it]][it];
3726       sumdy2 += dy[best[iter][it]][it]*dy[best[iter][it]][it];
3727       sumx   += dx;
3728       sumx2  += dx*dx;
3729       sumxy  += dx*dy[best[iter][it]][it];
3730       mpads  += cl[best[iter][it]][it]->GetNPads();
3731       if ((cl[best[iter][it]][it]->GetLabel(0) == label) ||
3732           (cl[best[iter][it]][it]->GetLabel(1) == label) ||
3733           (cl[best[iter][it]][it]->GetLabel(2) == label)) {
3734         ngood[iter]++;
3735       }
3736       else {
3737         nbad[iter]++;
3738       }
3739
3740     }
3741
3742     //
3743     // calculates line parameters
3744     //
3745     Double_t det  = sum*sumx2 - sumx*sumx;
3746     angle[iter]   = (sum*sumxy - sumx*sumdy) / det;
3747     mean[iter]    = (sumx2*sumdy - sumx*sumxy) / det;
3748     meanz[iter]   = sumz / sum;
3749     moffset[iter] = sumdy / sum;
3750     mpads        /= sum;   // Mean number of pads
3751
3752     Double_t sigma2 = 0.0; // Normalized residuals - for line fit
3753     Double_t sigma1 = 0.0; // Normalized residuals - constant fit
3754
3755     for (Int_t it = 0; it < t1 - t0; it++) {
3756       if (!cl[best[iter][it]][it]) {
3757         continue;
3758       }
3759       Double_t dx  = x[it] - xmean;
3760       Double_t ytr = mean[iter] + angle[iter] * dx;
3761       sigma2 += (dy[best[iter][it]][it] - ytr)
3762               * (dy[best[iter][it]][it] - ytr);
3763       sigma1 += (dy[best[iter][it]][it] - moffset[iter])
3764               * (dy[best[iter][it]][it] - moffset[iter]);
3765       sum++;
3766     }
3767     sigma2          /= (sum - 2);                  // Normalized residuals
3768     sigma1          /= (sum - 1);                  // Normalized residuals
3769     smean[iter]      = sigma2 * (sumx2 / det);     // Estimated error2 of mean
3770     sangle[iter]     = sigma2 * ( sum  / det);     // Estimated error2 of angle
3771     smeanangle[iter] = sigma2 * (-sumx / det);     // Correlation
3772     sigmas[iter]     = TMath::Sqrt(sigma1);
3773     smoffset[iter]   = (sigma1 / sum) + 0.01*0.01; // Sigma of mean offset + unisochronity sigma
3774
3775     //
3776     // Iterative choice of "better path"
3777     //
3778     for (Int_t it = 0; it < t1 - t0; it++) {
3779
3780       if (!cl[best[iter][it]][it]) {
3781         continue;
3782       }
3783
3784       // Add unisochronity + angular effect contribution
3785       Double_t sigmatr2 = smoffset[iter] + 0.5*tany*tany;
3786       Double_t sweight  = 1.0/sigmatr2 + 1.0/track->GetSigmaY2();
3787       Double_t weighty  = (moffset[iter] / sigmatr2) / sweight;             // Weighted mean
3788       Double_t sigmacl  = TMath::Sqrt(sigma1*sigma1 + track->GetSigmaY2());
3789       Double_t mindist  = 100000.0;
3790       Int_t    ihbest   = 0;
3791
3792       for (Int_t ih = 0; ih < 10; ih++) {
3793         if (!cl[ih][it]) {
3794           break;
3795         }
3796         Double_t dist2 = (dy[ih][it] - weighty) / sigmacl;
3797         dist2 *= dist2; // Chi2 distance
3798         if (dist2 < mindist) {
3799           mindist = dist2;
3800           ihbest  = ih;
3801         }
3802       }
3803       best[iter+1][it] = ihbest;
3804     }
3805
3806     //
3807     // Update best hypothesy if better chi2 according tracklet position and angle
3808     //
3809     Double_t sy2 = smean[iter]  + track->GetSigmaY2();
3810     Double_t sa2 = sangle[iter] + track->GetSigmaSnp2(); // track->fCee;
3811     Double_t say = track->GetSigmaSnpY();                // track->fCey;
3812     //Double_t chi20 = mean[bestiter]*mean[bestiter ] / sy2+angle[bestiter]*angle[bestiter]/sa2;
3813     //Double_t chi21 = mean[iter]*mean[iter]          / sy2+angle[iter]*angle[iter]/sa2;
3814
3815     Double_t detchi    = sy2*sa2 - say*say;
3816     Double_t invers[3] = {sa2/detchi,sy2/detchi,-say/detchi}; // Inverse value of covariance matrix
3817
3818     Double_t chi20 = mean[bestiter]  * mean[bestiter]       * invers[0]
3819                    + angle[bestiter] * angle[bestiter]      * invers[1]
3820                    + 2.0 * mean[bestiter] * angle[bestiter] * invers[2];
3821     Double_t chi21 = mean[iter]  * mean[iter]               * invers[0]
3822                    + angle[iter] * angle[iter]              * invers[1]
3823                    + 2.0 * mean[iter] * angle[iter]         * invers[2];
3824     tchi2s[iter]   = chi21;
3825
3826     if ((changes[iter] <= changes[bestiter]) &&
3827         (chi21 < chi20)) {
3828       bestiter = iter;
3829     }
3830
3831   }
3832
3833   //
3834   // Set clusters
3835   //
3836   Double_t sigma2     = sigmas[0]; // Choose as sigma from 0 iteration
3837   Short_t  maxpos     = -1;
3838   Float_t  maxcharge  =  0.0;
3839   Short_t  maxpos4    = -1;
3840   Float_t  maxcharge4 =  0.0;
3841   Short_t  maxpos5    = -1;
3842   Float_t  maxcharge5 =  0.0;
3843
3844   //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3845   //if (tchi2s[bestiter]>25.) sigma2=1000.;  // dont'accept
3846
3847   Double_t exB         = AliTRDcalibDB::Instance()->GetOmegaTau(
3848                           AliTRDcalibDB::Instance()->GetVdrift(0,0,0));
3849   Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
3850   if (mpads > 3.5) {
3851     expectederr += (mpads - 3.5) * 0.04;
3852   }
3853   if (changes[bestiter] > 1) {
3854     expectederr += changes[bestiter] * 0.01;
3855   }
3856   expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0;
3857   //if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
3858   //expectederr+=10000;
3859
3860   for (Int_t it = 0; it < t1 - t0; it++) {
3861
3862     if (!cl[best[bestiter][it]][it]) {
3863       continue;
3864     }
3865
3866     cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // Set cluster error
3867     if (!cl[best[bestiter][it]][it]->IsUsed()) {
3868       cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY());
3869       //cl[best[bestiter][it]][it]->Use();
3870     }
3871
3872     // Time bins with maximal charge
3873     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
3874       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3875       maxpos    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3876     }
3877
3878     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3879       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
3880         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3881         maxpos4    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3882       }
3883     }
3884
3885     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3886       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
3887         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3888         maxpos5    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3889       }
3890     }
3891
3892     // Time bins with maximal charge
3893     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge) {
3894       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3895       maxpos    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3896     }
3897
3898     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3899       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
3900         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3901         maxpos4    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3902       }
3903     }
3904
3905     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3906       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
3907         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3908         maxpos5    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3909       }
3910     }
3911
3912     clusters[it+t0] = indexes[best[bestiter][it]][it];
3913
3914     // Still needed ????
3915     //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 &&
3916     //cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0]
3917     // = indexes[best[bestiter][it]][it];    //Test
3918
3919   }
3920
3921   //
3922   // Set tracklet parameters
3923   //
3924   Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01;
3925   if (mpads > 3.5) {
3926     trackleterr2 += (mpads - 3.5) * 0.04;
3927   }
3928   trackleterr2 += changes[bestiter] * 0.01;
3929   trackleterr2 *= TMath::Max(14.0 - nfound,1.0);
3930   trackleterr2 += 0.2 * (tany-exB)*(tany-exB);
3931
3932   // Set tracklet parameters
3933   tracklet.Set(xmean
3934               ,track2.GetY() + moffset[bestiter]
3935               ,meanz[bestiter]
3936               ,track2.GetAlpha()
3937               ,trackleterr2);
3938   tracklet.SetTilt(h01);
3939   tracklet.SetP0(mean[bestiter]);
3940   tracklet.SetP1(angle[bestiter]);
3941   tracklet.SetN(nfound);
3942   tracklet.SetNCross(changes[bestiter]);
3943   tracklet.SetPlane(plane);
3944   tracklet.SetSigma2(expectederr);
3945   tracklet.SetChi2(tchi2s[bestiter]);
3946   tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
3947   track->fTracklets[plane] = tracklet;
3948   track->fNWrong += nbad[0];
3949
3950   //
3951   // Debuging part
3952   //
3953   TClonesArray array0("AliTRDcluster");
3954   TClonesArray array1("AliTRDcluster");
3955   array0.ExpandCreateFast(t1 - t0 + 1);
3956   array1.ExpandCreateFast(t1 - t0 + 1);
3957   TTreeSRedirector &cstream = *fDebugStreamer;
3958   AliTRDcluster dummy;
3959   Double_t dy0[100];
3960   Double_t dyb[100];
3961
3962   for (Int_t it = 0; it < t1 - t0; it++) {
3963     dy0[it] = dy[0][it];
3964     dyb[it] = dy[best[bestiter][it]][it];
3965     if (cl[0][it]) {
3966       new(array0[it]) AliTRDcluster(*cl[0][it]);
3967     }
3968     else {
3969       new(array0[it]) AliTRDcluster(dummy);
3970     }
3971     if(cl[best[bestiter][it]][it]) {
3972       new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
3973     }
3974     else{
3975       new(array1[it]) AliTRDcluster(dummy);
3976     }
3977   }
3978   TGraph graph0(t1-t0,x,dy0);
3979   TGraph graph1(t1-t0,x,dyb);
3980   TGraph graphy(t1-t0,x,yt);
3981   TGraph graphz(t1-t0,x,zt);
3982
3983   if (AliTRDReconstructor::StreamLevel() > 0) {
3984     cstream << "tracklet"
3985             << "track.="      << track              // Track parameters
3986             << "tany="        << tany               // Tangent of the local track angle
3987             << "xmean="       << xmean              // Xmean - reference x of tracklet
3988             << "tilt="        << h01                // Tilt angle
3989             << "nall="        << nall               // Number of foundable clusters
3990             << "nfound="      << nfound             // Number of found clusters
3991             << "clfound="     << clfound            // Total number of found clusters in road
3992             << "mpads="       << mpads              // Mean number of pads per cluster
3993             << "plane="       << plane              // Plane number
3994             << "detector="    << detector           // Detector number
3995             << "road="        << road               // The width of the used road
3996             << "graph0.="     << &graph0            // x - y = dy for closest cluster
3997             << "graph1.="     << &graph1            // x - y = dy for second closest cluster
3998             << "graphy.="     << &graphy            // y position of the track
3999             << "graphz.="     << &graphz            // z position of the track
4000             //<< "fCl.="        << &array0            // closest cluster
4001             //<< "fCl2.="       << &array1            // second closest cluster
4002             << "maxpos="      << maxpos             // Maximal charge postion
4003             << "maxcharge="   << maxcharge          // Maximal charge
4004             << "maxpos4="     << maxpos4            // Maximal charge postion - after bin 4
4005             << "maxcharge4="  << maxcharge4         // Maximal charge         - after bin 4
4006             << "maxpos5="     << maxpos5            // Maximal charge postion - after bin 5
4007             << "maxcharge5="  << maxcharge5         // Maximal charge         - after bin 5
4008             << "bestiter="    << bestiter           // Best iteration number
4009             << "tracklet.="   << &tracklet          // Corrspond to the best iteration
4010             << "tchi20="      << tchi2s[0]          // Chi2 of cluster in the 0 iteration
4011             << "tchi2b="      << tchi2s[bestiter]   // Chi2 of cluster in the best  iteration
4012             << "sigmas0="     << sigmas[0]          // Residuals sigma
4013             << "sigmasb="     << sigmas[bestiter]   // Residulas sigma
4014             << "ngood0="      << ngood[0]           // Number of good clusters in 0 iteration
4015             << "nbad0="       << nbad[0]            // Number of bad clusters in 0 iteration
4016             << "ngoodb="      << ngood[bestiter]    //                        in best iteration
4017             << "nbadb="       << nbad[bestiter]     //                        in best iteration
4018             << "changes0="    << changes[0]         // Changes of pardrows in iteration number 0
4019             << "changesb="    << changes[bestiter]  // Changes of pardrows in best iteration
4020             << "moffset0="    << moffset[0]         // Offset fixing angle in iter=0
4021             << "smoffset0="   << smoffset[0]        // Sigma of offset fixing angle in iter=0
4022             << "moffsetb="    << moffset[bestiter]  // Offset fixing angle in iter=best
4023             << "smoffsetb="   << smoffset[bestiter] // Sigma of offset fixing angle in iter=best
4024             << "mean0="       << mean[0]            // Mean dy in iter=0;
4025             << "smean0="      << smean[0]           // Sigma of mean dy in iter=0
4026             << "meanb="       << mean[bestiter]     // Mean dy in iter=best
4027             << "smeanb="      << smean[bestiter]    // Sigma of mean dy in iter=best
4028             << "angle0="      << angle[0]           // Angle deviation in the iteration number 0
4029             << "sangle0="     << sangle[0]          // Sigma of angular deviation in iteration number 0
4030             << "angleb="      << angle[bestiter]    // Angle deviation in the best iteration
4031             << "sangleb="     << sangle[bestiter]   // Sigma of angle deviation in the best iteration
4032             << "expectederr=" << expectederr        // Expected error of cluster position
4033             << "\n";
4034   }
4035
4036   return nfound;
4037
4038 }
4039
4040 //_____________________________________________________________________________
4041 Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
4042                         , Int_t *outlist, Bool_t down)
4043 {
4044   //
4045   // Sort eleements according occurancy
4046   // The size of output array has is 2*n
4047   //
4048
4049   Int_t *sindexS = new Int_t[n];   // Temporary array for sorting
4050   Int_t *sindexF = new Int_t[2*n];
4051   for (Int_t i = 0; i < n; i++) {
4052     sindexF[i] = 0;
4053   }
4054
4055   TMath::Sort(n,inlist,sindexS,down);
4056
4057   Int_t last     = inlist[sindexS[0]];
4058   Int_t val      = last;
4059   sindexF[0]     = 1;
4060   sindexF[0+n]   = last;
4061   Int_t countPos = 0;
4062
4063   // Find frequency
4064   for (Int_t i = 1; i < n; i++) {
4065     val = inlist[sindexS[i]];
4066     if (last == val) {
4067       sindexF[countPos]++;
4068     }
4069     else {
4070       countPos++;
4071       sindexF[countPos+n] = val;
4072       sindexF[countPos]++;
4073       last                = val;
4074     }
4075   }
4076   if (last == val) {
4077     countPos++;
4078   }
4079
4080   // Sort according frequency
4081   TMath::Sort(countPos,sindexF,sindexS,kTRUE);
4082
4083   for (Int_t i = 0; i < countPos; i++) {
4084     outlist[2*i  ] = sindexF[sindexS[i]+n];
4085     outlist[2*i+1] = sindexF[sindexS[i]];
4086   }
4087
4088   delete [] sindexS;
4089   delete [] sindexF;
4090
4091   return countPos;
4092
4093 }
4094
4095 //_____________________________________________________________________________
4096 AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
4097 {
4098   //
4099   // Register a seed
4100   //
4101
4102   Double_t alpha = AliTRDgeometry::GetAlpha();
4103   Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
4104   Double_t c[15];
4105
4106   c[ 0] = 0.2;
4107   c[ 1] = 0.0; c[ 2] = 2.0;
4108   c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
4109   c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0;  c[ 9] = 0.1;
4110   c[10] = 0.0; c[11] = 0.0; c[12] = 0.0;  c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
4111
4112   Int_t index = 0;
4113   AliTRDcluster *cl = 0;
4114
4115   for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
4116     if (seeds[ilayer].IsOK()) {
4117       for (Int_t itime = 22; itime > 0; itime--) {
4118         if (seeds[ilayer].GetIndexes(itime) > 0) {
4119           index = seeds[ilayer].GetIndexes(itime);
4120           cl    = seeds[ilayer].GetClusters(itime);
4121           break;
4122         }
4123       }
4124     }
4125     if (index > 0) {
4126       break;
4127     }
4128   }
4129   if (cl == 0) {
4130     return 0;
4131   }
4132
4133   AliTRDtrack *track = new AliTRDtrack(cl
4134                                       ,index
4135                                       ,&params[1]
4136                                       ,c
4137                                       ,params[0]
4138                                       ,params[6]*alpha+shift);
4139   track->PropagateTo(params[0]-5.0);
4140   track->ResetCovariance(1);
4141
4142   Int_t rc = FollowBackProlongation(*track);
4143   if (rc < 30) {
4144     delete track;
4145     track = 0;
4146   }
4147   else {
4148     track->CookdEdx();
4149     CookdEdxTimBin(*track);
4150     CookLabel(track,0.9);
4151   }
4152
4153   return track;
4154
4155 }