TRD/AliTRDtracker.cxx

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