TRD/AliTRDtracker.cxx

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