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