HLT/TRD/AliTRDtrackerHLT.cxx

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