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