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