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