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