TRD/AliTRDtracker.cxx

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