TRD/AliTRDtracker.cxx

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