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