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