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