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