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