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