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