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