TRD/AliTRDtracker.cxx

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