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