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