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