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     AliErrorGeneral("AliTRDtrackingSector::Ctor"
2748                    ,"Could not get common parameters\n");
2749     return;
2750   }
2751
2752   for (Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2753
2754     ymax          = fGeom->GetChamberWidth(plane) / 2.0;
2755     padPlane      = commonParam->GetPadPlane(plane,0);
2756     ymaxsensitive = (padPlane->GetColSize(1) * padPlane->GetNcols() - 4.0) / 2.0;
2757
2758     for (Int_t ch = 0; ch < kNchambers; ch++) {
2759       zmax[ch]          = fGeom->GetChamberLength(plane,ch) / 2.0;
2760       Float_t pad       = padPlane->GetRowSize(1);
2761       Float_t row0      = commonParam->GetRow0(plane,ch,0);
2762       Int_t   nPads     = commonParam->GetRowMax(plane,ch,0);
2763       zmaxsensitive[ch] = Float_t(nPads) * pad / 2.0;
2764       zc[ch]            = -(pad * nPads) / 2.0 + row0;
2765     }
2766
2767     dx        = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
2768               / AliTRDcalibDB::Instance()->GetSamplingFrequency();
2769     rho       = 0.00295 * 0.85; //????
2770     radLength = 11.0;
2771
2772     Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
2773     //Double_t xbottom = x0 - dxDrift;
2774     //Double_t xtop    = x0 + dxAmp;
2775
2776     Int_t nTimeBins =  AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
2777     for (Int_t iTime = 0; iTime < nTimeBins; iTime++) {
2778
2779       Double_t xlayer = iTime * dx - dxAmp;
2780       //if (xlayer<0) xlayer = dxAmp / 2.0;
2781       x = x0 - xlayer;
2782
2783       tbIndex = CookTimeBinIndex(plane,iTime);
2784       ppl     = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex,plane);
2785       ppl->SetYmax(ymax,ymaxsensitive);
2786       ppl->SetZ(zc,zmax,zmaxsensitive);
2787       ppl->SetHoles(holes);
2788       InsertLayer(ppl);
2789
2790     }
2791
2792   }
2793
2794   MapTimeBinLayers();
2795
2796   delete [] zc;
2797   delete [] zmax;
2798   delete [] zmaxsensitive;
2799
2800 }
2801
2802 //_____________________________________________________________________________
2803 AliTRDtracker::AliTRDtrackingSector
2804              ::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
2805   :fN(0)
2806   ,fGeom(0)
2807   ,fGeomSector(0)
2808 {
2809   //
2810   // Copy constructor
2811   //
2812
2813 }
2814
2815 //_____________________________________________________________________________
2816 Int_t  AliTRDtracker::AliTRDtrackingSector
2817                     ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
2818 {
2819   //
2820   // Depending on the digitization parameters calculates "global"
2821   // time bin index for timebin <localTB> in plane <plane>
2822   //
2823   //
2824
2825   Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
2826   Int_t gtb        = (plane + 1) * tbPerPlane - localTB - 1;
2827
2828   if (localTB < 0) return -1;
2829   if (gtb     < 0) return -1;
2830
2831   return gtb;
2832
2833 }
2834
2835 //_____________________________________________________________________________
2836 void AliTRDtracker::AliTRDtrackingSector
2837                   ::MapTimeBinLayers()
2838 {
2839   //
2840   // For all sensitive time bins sets corresponding layer index
2841   // in the array fTimeBins
2842   //
2843
2844   Int_t index;
2845
2846   for (Int_t i = 0; i < fN; i++) {
2847
2848     index = fLayers[i]->GetTimeBinIndex();
2849
2850     if(index < 0) continue;
2851     if(index >= (Int_t) kMaxTimeBinIndex) {
2852       AliWarningGeneral("AliTRDtrackingSector::MapTimeBinLayers()"
2853                        ,Form("Index %d exceeds allowed maximum of %d!\n"
2854                        ,index,kMaxTimeBinIndex-1));
2855       continue;
2856     }
2857
2858     fTimeBinIndex[index] = i;
2859
2860   }
2861
2862 }
2863
2864 //_____________________________________________________________________________
2865 Int_t AliTRDtracker::AliTRDtrackingSector
2866                    ::GetLayerNumber(Double_t x) const
2867 {
2868   //
2869   // Returns the number of time bin which in radial position is closest to <x>
2870   //
2871
2872   if (x >= fLayers[fN-1]->GetX()) return fN - 1;
2873   if (x <= fLayers[0]->GetX())    return 0;
2874
2875   Int_t b = 0;
2876   Int_t e = fN - 1;
2877   Int_t m = (b + e) / 2;
2878
2879   for (; b < e; m = (b + e)/ 2) {
2880     if (x > fLayers[m]->GetX()) {
2881       b = m + 1;
2882     }
2883     else {
2884       e = m;
2885     }
2886   }
2887
2888   if (TMath::Abs(x - fLayers[m]->GetX()) >
2889       TMath::Abs(x - fLayers[m+1]->GetX())) {
2890     return m+1;
2891   }
2892   else {
2893     return m;
2894   }
2895
2896 }
2897
2898 //_____________________________________________________________________________
2899 Int_t AliTRDtracker::AliTRDtrackingSector
2900                    ::GetInnerTimeBin() const
2901 {
2902   //
2903   // Returns number of the innermost SENSITIVE propagation layer
2904   //
2905
2906   return GetLayerNumber(0);
2907
2908 }
2909
2910 //_____________________________________________________________________________
2911 Int_t AliTRDtracker::AliTRDtrackingSector
2912                    ::GetOuterTimeBin() const
2913 {
2914   //
2915   // Returns number of the outermost SENSITIVE time bin
2916   //
2917
2918   return GetLayerNumber(GetNumberOfTimeBins() - 1);
2919
2920 }
2921
2922 //_____________________________________________________________________________
2923 Int_t AliTRDtracker::AliTRDtrackingSector
2924                    ::GetNumberOfTimeBins() const
2925 {
2926   //
2927   // Returns number of SENSITIVE time bins
2928   //
2929
2930   Int_t tb;
2931   Int_t layer;
2932
2933   for (tb = kMaxTimeBinIndex - 1; tb >= 0; tb--) {
2934     layer = GetLayerNumber(tb);
2935     if (layer >= 0) break;
2936   }
2937
2938   return tb + 1;
2939
2940 }
2941
2942 //_____________________________________________________________________________
2943 void AliTRDtracker::AliTRDtrackingSector
2944                   ::InsertLayer(AliTRDpropagationLayer* pl)
2945 {
2946   //
2947   // Insert layer <pl> in fLayers array.
2948   // Layers are sorted according to X coordinate.
2949   //
2950
2951   if (fN == ((Int_t) kMaxLayersPerSector)) {
2952     AliWarningGeneral("AliTRDtrackingSector::InsertLayer"
2953                      ,"Too many layers !\n");
2954     return;
2955   }
2956
2957   if (fN == 0) {
2958     fLayers[fN++] = pl;
2959     return;
2960   }
2961
2962   Int_t i = Find(pl->GetX());
2963   memmove(fLayers + i + 1,fLayers + i
2964          ,(fN-i)*sizeof(AliTRDpropagationLayer*));
2965   fLayers[i] = pl;
2966   fN++;
2967
2968 }
2969
2970 //_____________________________________________________________________________
2971 Int_t AliTRDtracker::AliTRDtrackingSector
2972                    ::Find(Double_t x) const
2973 {
2974   //
2975   // Returns index of the propagation layer nearest to X
2976   //
2977
2978   if (x <= fLayers[0]->GetX())    return 0;
2979   if (x >  fLayers[fN-1]->GetX()) return fN;
2980
2981   Int_t b = 0;
2982   Int_t e = fN - 1;
2983   Int_t m = (b + e) / 2;
2984
2985   for (; b < e; m = (b + e) / 2) {
2986     if (x > fLayers[m]->GetX()) {
2987       b = m + 1;
2988     }
2989     else {
2990       e = m;
2991     }
2992   }
2993
2994   return m;
2995
2996 }
2997
2998 //_____________________________________________________________________________
2999 void AliTRDtracker::AliTRDpropagationLayer
3000                   ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
3001 {
3002   //
3003   // Set centers and the width of sectors
3004   //
3005
3006   for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3007
3008     fZc[icham]            = center[icham];
3009     fZmax[icham]          = w[icham];
3010     fZmaxSensitive[icham] = wsensitive[icham];
3011
3012   }
3013
3014 }
3015
3016 //_____________________________________________________________________________
3017 void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3018 {
3019   //
3020   // Set centers and the width of sectors
3021   //
3022
3023   fHole = kFALSE;
3024
3025   for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
3026     fIsHole[icham] = holes[icham];
3027     if (holes[icham]) {
3028       fHole = kTRUE;
3029     }
3030   }
3031
3032 }
3033
3034 //_____________________________________________________________________________
3035 void AliTRDtracker::AliTRDpropagationLayer
3036                   ::InsertCluster(AliTRDcluster *c, UInt_t index)
3037 {
3038   //
3039   // Insert cluster in cluster array.
3040   // Clusters are sorted according to Y coordinate.
3041   //
3042
3043   if (fTimeBinIndex < 0) {
3044     AliErrorGeneral("AliTRDpropagationLayer::InsertCluster"
3045                    ,"Attempt to insert cluster into non-sensitive time bin!\n");
3046     return;
3047   }
3048
3049   if (fN == (Int_t) kMaxClusterPerTimeBin) {
3050     AliErrorGeneral("AliTRDpropagationLayer::InsertCluster"
3051                    ,"Too many clusters !\n");
3052     return;
3053   }
3054
3055   if (fN == 0) {
3056     fIndex[0]       = index;
3057     fClusters[fN++] = c;
3058     return;
3059   }
3060
3061   Int_t i = Find(c->GetY());
3062   memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3063   memmove(fIndex   +i+1,fIndex   +i,(fN-i)*sizeof(UInt_t));
3064   fIndex[i]    = index;
3065   fClusters[i] = c;
3066   fN++;
3067
3068 }
3069
3070 //_____________________________________________________________________________
3071 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
3072 {
3073   //
3074   // Returns index of the cluster nearest in Y
3075   //
3076
3077   if (fN <= 0)                       return 0;
3078   if (y  <= fClusters[0]->GetY())    return 0;
3079   if (y  >  fClusters[fN-1]->GetY()) return fN;
3080
3081   Int_t b = 0;
3082   Int_t e = fN - 1;
3083   Int_t m = (b + e) / 2;
3084
3085   for (; b < e; m = (b + e) / 2) {
3086     if (y > fClusters[m]->GetY()) {
3087       b = m + 1;
3088     }
3089     else {
3090       e = m;
3091     }
3092   }
3093
3094   return m;
3095
3096 }
3097
3098 //_____________________________________________________________________________
3099 Int_t AliTRDtracker::AliTRDpropagationLayer
3100                    ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
3101                                       , Float_t maxroadz) const
3102 {
3103   //
3104   // Returns index of the cluster nearest to the given y,z
3105   //
3106
3107   Int_t   index   = -1;
3108   Int_t   maxn    = fN;
3109   Float_t mindist = maxroad;
3110
3111   for (Int_t i = Find(y - maxroad); i < maxn; i++) {
3112     AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
3113     Float_t ycl = c->GetY();
3114     if (ycl > y + maxroad) break;
3115     if (TMath::Abs(c->GetZ() - z) > maxroadz) continue;
3116     if (TMath::Abs(ycl - y)       < mindist) {
3117       mindist = TMath::Abs(ycl - y);
3118       index   = fIndex[i];
3119     }
3120   }
3121
3122   return index;
3123
3124 }
3125
3126 //_____________________________________________________________________________
3127 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c)
3128 {
3129   //
3130   // Returns correction factor for tilted pads geometry
3131   //
3132
3133   Int_t det   = c->GetDetector();
3134   Int_t plane = fGeom->GetPlane(det);
3135
3136   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
3137
3138   Double_t h01 = TMath::Tan(-TMath::Pi()/180.0 * padPlane->GetTiltingAngle());
3139
3140   if (fNoTilt) h01 = 0.0;
3141
3142   return h01;
3143
3144 }
3145
3146 //_____________________________________________________________________________
3147 void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
3148 {
3149   //
3150   // This is setting fdEdxPlane and fTimBinPlane
3151   // Sums up the charge in each plane for track TRDtrack and also get the
3152   // Time bin for Max. Cluster
3153   // Prashant Shukla (shukla@physi.uni-heidelberg.de)
3154   //
3155
3156   Double_t  clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
3157   Double_t  maxclscharge[AliESDtrack::kNPlane];
3158   Int_t     nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
3159   Int_t     timebin[AliESDtrack::kNPlane];
3160
3161   // Initialization of cluster charge per plane.
3162   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
3163     for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
3164       clscharge[iPlane][iSlice] = 0.0;
3165       nCluster[iPlane][iSlice]  = 0;
3166     }
3167   }
3168
3169   // Initialization of cluster charge per plane.
3170   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
3171     timebin[iPlane]      = -1;
3172     maxclscharge[iPlane] = 0.0;
3173   }
3174
3175   // Loop through all clusters associated to track TRDtrack
3176   Int_t nClus = TRDtrack.GetNumberOfClusters();  // from Kalmantrack
3177   for (Int_t iClus = 0; iClus < nClus; iClus++) {
3178     Double_t charge = TRDtrack.GetClusterdQdl(iClus);
3179     Int_t    index  = TRDtrack.GetClusterIndex(iClus);
3180     AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index);
3181     if (!pTRDcluster) continue;
3182     Int_t tb = pTRDcluster->GetLocalTimeBin();
3183     if (!tb)          continue;
3184     Int_t detector = pTRDcluster->GetDetector();
3185     Int_t iPlane   = fGeom->GetPlane(detector);
3186     Int_t iSlice   = tb * AliESDtrack::kNSlice / AliTRDtrack::kNtimeBins;
3187     clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge;
3188     if(charge > maxclscharge[iPlane]) {
3189       maxclscharge[iPlane] = charge;
3190       timebin[iPlane]      = tb;
3191     }
3192     nCluster[iPlane][iSlice]++;
3193   }
3194
3195   // Setting the fdEdxPlane and fTimBinPlane variabales
3196   Double_t totalCharge = 0;
3197   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
3198     for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
3199       if (nCluster[iPlane][iSlice]) {
3200         clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
3201       }
3202       TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice);
3203       totalCharge = totalCharge + clscharge[iPlane][iSlice];
3204     }
3205     TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);
3206   }
3207
3208   // Still needed ????
3209   //  Int_t i;
3210   //  Int_t nc=TRDtrack.GetNumberOfClusters();
3211   //  Float_t dedx=0;
3212   //  for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
3213   //  dedx /= nc;
3214   //  for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3215   //    TRDtrack.SetPIDsignals(dedx, iPlane);
3216   //    TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
3217   //  }
3218
3219 }
3220
3221 //_____________________________________________________________________________
3222 Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
3223                                 , AliTRDtrack *track
3224                                 , Int_t *clusters
3225                                 , AliTRDtracklet &tracklet)
3226 {
3227   //
3228   // Try to find nearest clusters to the track in timebins from t0 to t1
3229   //
3230   // correction coefficients
3231   // - depend on TRD parameters
3232   // - to be changed according it
3233   //
3234
3235   Double_t x[100];
3236   Double_t yt[100];
3237   Double_t zt[100];
3238   Double_t xmean = 0.0;       // Reference x
3239   Double_t dz[10][100];
3240   Double_t dy[10][100];
3241   Float_t  zmean[100];
3242   Float_t  nmean[100];
3243   Int_t    clfound = 0;
3244   Int_t    indexes[10][100];  // Indexes of the clusters in the road
3245   Int_t    best[10][100];     // Index of best matching cluster
3246
3247   AliTRDcluster *cl[10][100]; // Pointers to the clusters in the road
3248
3249   for (Int_t it = 0; it < 100; it++) {
3250
3251     x[it]        =  0.0;
3252     yt[it]       =  0.0;
3253     zt[it]       =  0.0;
3254     clusters[it] = -2;
3255     zmean[it]    =  0.0;
3256     nmean[it]    =  0.0;
3257
3258     for (Int_t ih = 0; ih < 10; ih++) {
3259       indexes[ih][it] = -2;              // Reset indexes1
3260       cl[ih][it]      =  0;
3261       dz[ih][it]      = -100.0;
3262       dy[ih][it]      = -100.0;
3263       best[ih][it]    =  0;
3264     }
3265
3266   }
3267
3268   AliTRDtrack track2(*track);
3269   Double_t x0        = track->GetX();
3270   Double_t sigmaz    = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
3271   Int_t    nall      = 0;
3272   Int_t    nfound    = 0;
3273   Double_t h01       = 0.0;
3274   Int_t    plane     = -1;
3275   Int_t    detector  = -1;
3276   Float_t  padlength = 0.0;
3277   Float_t  snpy      = track->GetSnp();
3278   Float_t  tany      = TMath::Sqrt(snpy*snpy / (1.0 - snpy*snpy));
3279
3280   if (snpy < 0.0) {
3281     tany *= -1.0;
3282   }
3283
3284   Double_t sy2       = ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
3285   Double_t sz2       = ExpectedSigmaZ2(x0,track->GetTgl());
3286   Double_t road      = 15.0 * TMath::Sqrt(track->GetSigmaY2() + sy2);
3287   if (road > 6.0) {
3288     road = 6.0;
3289   }
3290
3291   for (Int_t it = 0; it < t1-t0; it++) {
3292
3293     Double_t maxChi2[2] = { fgkMaxChi2, fgkMaxChi2 };
3294     AliTRDpropagationLayer &timeBin= *(fTrSec[sector]->GetLayer(it+t0));
3295     if (timeBin == 0) continue;  // No indexes1
3296
3297     Int_t maxn = timeBin;
3298     x[it]      = timeBin.GetX();
3299     track2.PropagateTo(x[it]);
3300     yt[it]     = track2.GetY();
3301     zt[it]     = track2.GetZ();
3302
3303     Double_t y    = yt[it];
3304     Double_t z    = zt[it];
3305     Double_t chi2 = 1000000.0;
3306     nall++;
3307
3308     //
3309     // Find 2 nearest cluster at given time bin
3310     //
3311     for (Int_t i = timeBin.Find(y-road); i < maxn; i++) {
3312
3313       AliTRDcluster *c= (AliTRDcluster *) (timeBin[i]);
3314       h01 = GetTiltFactor(c);
3315       if (plane < 0) {
3316         Int_t det = c->GetDetector();
3317         plane     = fGeom->GetPlane(det);
3318         padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
3319       }
3320
3321       if (c->GetY() > y+road) break;
3322       if (((c->GetZ() - z) * (c->GetZ() - z)) > (12.0 * sz2)) continue;
3323
3324       Double_t dist = TMath::Abs(c->GetZ() - z);
3325       if (dist > (0.5 * padlength + 6.0 * sigmaz)) continue;   // 6 sigma boundary cut
3326       Double_t cost = 0.0;
3327       // Sigma boundary cost function
3328       if (dist > (0.5 * padlength - sigmaz)) {
3329         cost = (dist - 0.5 * padlength) / (2.0 * sigmaz);
3330         if (cost > -1.0) {
3331           cost = (cost + 1.0) * (cost + 1.0);
3332         }
3333         else {
3334           cost = 0.0;
3335         }
3336       }
3337       // Still needed ????
3338       //Int_t label = TMath::Abs(track->GetLabel());
3339       //if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
3340       chi2 = track2.GetPredictedChi2(c,h01) + cost;
3341
3342       clfound++;
3343       if (chi2 > maxChi2[1]) continue;
3344       detector = c->GetDetector();
3345
3346       // Store the clusters in the road
3347       for (Int_t ih = 2; ih < 9; ih++) {
3348         if (cl[ih][it] == 0) {
3349           cl[ih][it]      = c;
3350           indexes[ih][it] = timeBin.GetIndex(i); // Index - 9 - reserved for outliers
3351           break;
3352         }
3353       }
3354
3355       if (chi2 < maxChi2[0]) {
3356         maxChi2[1]     = maxChi2[0];
3357         maxChi2[0]     = chi2;
3358         indexes[1][it] = indexes[0][it];
3359         cl[1][it]      = cl[0][it];
3360         indexes[0][it] = timeBin.GetIndex(i);
3361         cl[0][it]      = c;
3362         continue;
3363       }
3364       maxChi2[1]     = chi2;
3365       cl[1][it]      = c;
3366       indexes[1][it] = timeBin.GetIndex(i);
3367
3368     }
3369
3370     if (cl[0][it]) {
3371       nfound++;
3372       xmean += x[it];
3373     }
3374
3375   }
3376
3377   if (nfound < 4) return 0;
3378   xmean /= Float_t(nfound);    // Middle x
3379   track2.PropagateTo(xmean);   // Propagate track to the center
3380
3381   //
3382   // Choose one of the variants
3383   //
3384   Int_t    changes[10];
3385   Float_t  sumz   = 0;
3386   Float_t  sum    = 0;
3387   Double_t sumdy  = 0;
3388   Double_t sumdy2 = 0;
3389   Double_t sumx   = 0;
3390   Double_t sumxy  = 0;
3391   Double_t sumx2  = 0;
3392   Double_t mpads  = 0;
3393
3394   Int_t    ngood[10];
3395   Int_t    nbad[10];
3396
3397   Double_t meanz[10];
3398   Double_t moffset[10];    // Mean offset
3399   Double_t mean[10];       // Mean value
3400   Double_t angle[10];      // Angle
3401
3402   Double_t smoffset[10];   // Sigma of mean offset
3403   Double_t smean[10];      // Sigma of mean value
3404   Double_t sangle[10];     // Sigma of angle
3405   Double_t smeanangle[10]; // Correlation
3406
3407   Double_t sigmas[10];
3408   Double_t tchi2s[10];     // Chi2s for tracklet
3409
3410   for (Int_t it = 0; it < 10; it++) {
3411
3412     ngood[it]      = 0;
3413     nbad[it]       = 0;
3414
3415     meanz[it]      = 0.0;
3416     moffset[it]    = 0.0;    // Mean offset
3417     mean[it]       = 0.0;    // Mean value
3418     angle[it]      = 0.0;    // Angle
3419
3420     smoffset[it]   = 1.0e10; // Sigma of mean offset
3421     smean[it]      = 1.0e10; // Sigma of mean value
3422     sangle[it]     = 1.0e10; // Sigma of angle
3423     smeanangle[it] = 0.0;    // Correlation
3424
3425     sigmas[it]     = 1.0e10;
3426     tchi2s[it]     = 1.0e10; // Chi2s for tracklet
3427
3428   }
3429
3430   //
3431   // Calculate zmean
3432   //
3433   for (Int_t it = 0; it < t1-t0; it++) {
3434
3435     if (!cl[0][it]) continue;
3436
3437     for (Int_t dt = -3; dt <= 3; dt++) {
3438       if (it+dt <     0) continue;
3439       if (it+dt > t1-t0) continue;
3440       if (!cl[0][it+dt]) continue;
3441       zmean[it] += cl[0][it+dt]->GetZ();
3442       nmean[it] += 1.0;
3443     }
3444     zmean[it] /= nmean[it];
3445
3446   }
3447
3448   for (Int_t it = 0; it < t1-t0; it++) {
3449
3450     best[0][it] = 0;
3451     for (Int_t ih = 0; ih < 10; ih++) {
3452
3453       dz[ih][it] = -100.0;
3454       dy[ih][it] = -100.0;
3455       if (!cl[ih][it]) continue;
3456       Double_t xcluster = cl[ih][it]->GetX();
3457       Double_t ytrack;
3458       Double_t ztrack;
3459       track2.GetProlongation(xcluster,ytrack,ztrack);
3460       dz[ih][it] = cl[ih][it]->GetZ() - ztrack;                   // Calculate z-distance from track
3461       dy[ih][it] = cl[ih][it]->GetY() + dz[ih][it]*h01  - ytrack; // Calculate y-distance from track
3462
3463     }
3464
3465     // Minimize changes
3466     if (!cl[0][it]) continue;
3467     if ((TMath::Abs(cl[0][it]->GetZ() - zmean[it]) > padlength * 0.8) &&
3468         (cl[1][it])) {
3469       if (TMath::Abs(cl[1][it]->GetZ() - zmean[it]) < padlength * 0.5) {
3470         best[0][it] = 1;
3471       }
3472     }
3473   }
3474
3475   //
3476   // Iterative choice of "best path"
3477   //
3478   Int_t label    = TMath::Abs(track->GetLabel());
3479   Int_t bestiter = 0;
3480
3481   for (Int_t iter = 0; iter < 9; iter++) {
3482
3483     changes[iter] = 0;
3484     sumz          = 0;
3485     sum           = 0;
3486     sumdy         = 0;
3487     sumdy2        = 0;
3488     sumx          = 0;
3489     sumx2         = 0;
3490     sumxy         = 0;
3491     mpads         = 0;
3492     ngood[iter]   = 0;
3493     nbad[iter]    = 0;
3494
3495     // Linear fit
3496     for (Int_t it = 0; it < t1-t0; it++) {
3497
3498       if (!cl[best[iter][it]][it]) continue;
3499
3500       // Calculates pad-row changes
3501       Double_t zbefore = cl[best[iter][it]][it]->GetZ();
3502       Double_t zafter  = cl[best[iter][it]][it]->GetZ();
3503       for (Int_t itd = it-1; itd >= 0; itd--) {
3504         if (cl[best[iter][itd]][itd]) {
3505           zbefore = cl[best[iter][itd]][itd]->GetZ();
3506           break;
3507         }
3508       }
3509       for (Int_t itd = it+1; itd < t1-t0; itd++) {
3510         if (cl[best[iter][itd]][itd]) {
3511           zafter  = cl[best[iter][itd]][itd]->GetZ();
3512           break;
3513         }
3514       }
3515       if ((TMath::Abs(cl[best[iter][it]][it]->GetZ() - zbefore) > 0.1) &&
3516           (TMath::Abs(cl[best[iter][it]][it]->GetZ() -  zafter) > 0.1)) {
3517         changes[iter]++;
3518       }
3519
3520       // Distance to reference x
3521       Double_t dx = x[it]-xmean;
3522       sumz   += cl[best[iter][it]][it]->GetZ();
3523       sum++;
3524       sumdy  += dy[best[iter][it]][it];
3525       sumdy2 += dy[best[iter][it]][it]*dy[best[iter][it]][it];
3526       sumx   += dx;
3527       sumx2  += dx*dx;
3528       sumxy  += dx*dy[best[iter][it]][it];
3529       mpads  += cl[best[iter][it]][it]->GetNPads();
3530       if ((cl[best[iter][it]][it]->GetLabel(0) == label) ||
3531           (cl[best[iter][it]][it]->GetLabel(1) == label) ||
3532           (cl[best[iter][it]][it]->GetLabel(2) == label)) {
3533         ngood[iter]++;
3534       }
3535       else {
3536         nbad[iter]++;
3537       }
3538
3539     }
3540
3541     //
3542     // Calculates line parameters
3543     //
3544     Double_t det  = sum * sumx2 - sumx * sumx;
3545     angle[iter]   = (sum * sumxy - sumx * sumdy) / det;
3546     mean[iter]    = (sumx2 * sumdy - sumx * sumxy) / det;
3547     meanz[iter]   = sumz  / sum;
3548     moffset[iter] = sumdy / sum;
3549     mpads        /= sum;   // Mean number of pads
3550
3551     Double_t sigma2 = 0.0; // Normalized residuals - for line fit
3552     Double_t sigma1 = 0.0; // Normalized residuals - constant fit
3553
3554     for (Int_t it = 0; it < t1-t0; it++) {
3555       if (!cl[best[iter][it]][it]) continue;
3556       Double_t dx  = x[it] - xmean;
3557       Double_t ytr = mean[iter] + angle[iter] * dx;
3558       sigma2 += (dy[best[iter][it]][it] - ytr)
3559               * (dy[best[iter][it]][it] - ytr);
3560       sigma1 += (dy[best[iter][it]][it] - moffset[iter])
3561               * (dy[best[iter][it]][it] - moffset[iter]);
3562       sum++;
3563     }
3564
3565     sigma2           /= (sum - 2.0);                // Normalized residuals
3566     sigma1           /= (sum - 1.0);                // Normalized residuals
3567
3568     smean[iter]       = sigma2 * (sumx2/det);       // Estimated error2 of mean
3569     sangle[iter]      = sigma2 * (  sum/det);       // Estimated error2 of angle
3570     smeanangle[iter]  = sigma2 * (-sumx/det);       // Correlation
3571
3572     sigmas[iter]      = TMath::Sqrt(sigma1);
3573     smoffset[iter]    = (sigma1 / sum) + 0.01*0.01; // Sigma of mean offset + unisochronity sigma
3574
3575     //
3576     // Iterative Choice of "better path"
3577     //
3578     for (Int_t it = 0; it < t1-t0; it++) {
3579
3580       if (!cl[best[iter][it]][it]) continue;
3581
3582       // Add unisochronity + angular effect contribution
3583       Double_t sigmatr2 = smoffset[iter] + 0.5*tany*tany;
3584       Double_t sweight  = 1.0 / sigmatr2 + 1.0 / track->GetSigmaY2();
3585       Double_t weighty  = (moffset[iter] / sigmatr2) / sweight; // Weighted mean
3586       Double_t sigmacl  = TMath::Sqrt(sigma1*sigma1 + track->GetSigmaY2());
3587       Double_t mindist  = 100000.0;
3588       Int_t    ihbest   = 0;
3589       for (Int_t ih = 0; ih < 10; ih++) {
3590         if (!cl[ih][it]) break;
3591         Double_t dist2  = (dy[ih][it] - weighty) / sigmacl;
3592         dist2          *= dist2; // Chi2 distance
3593         if (dist2 < mindist) {
3594           mindist = dist2;
3595           ihbest  = ih;
3596         }
3597       }
3598       best[iter+1][it] = ihbest;
3599
3600     }
3601
3602     //
3603     // Update best hypothesy if better chi2 according tracklet position and angle
3604     //
3605     Double_t sy2 = smean[iter]  + track->GetSigmaY2();
3606     Double_t sa2 = sangle[iter] + track->fCee;
3607     Double_t say = track->fCey;
3608     //Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2;
3609     //Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2;
3610
3611     Double_t detchi    = sy2*sa2 - say*say;
3612     // Inverse value of covariance matrix
3613     Double_t invers[3] = {  sa2/detchi,  sy2/detchi, -say/detchi};
3614
3615     Double_t chi20     =       mean[bestiter]*mean[bestiter]   * invers[0]
3616                        +       angle[bestiter]*angle[bestiter] * invers[1]
3617                        + 2.0 * mean[bestiter]*angle[bestiter]  * invers[2];
3618     Double_t chi21     =       mean[iter]*mean[iter]   * invers[0]
3619                        +       angle[iter]*angle[iter] * invers[1]
3620                        + 2.0 * mean[iter]*angle[iter]  * invers[2];
3621     tchi2s[iter]       = chi21;
3622
3623     if ((changes[iter] <= changes[bestiter]) &&
3624         (chi21         <  chi20)) {
3625       bestiter = iter;
3626     }
3627
3628   }
3629
3630   //
3631   // Set clusters
3632   //
3633   Double_t sigma2     = sigmas[0]; // Choose as sigma  from 0 iteration
3634   Short_t  maxpos     = -1;
3635   Float_t  maxcharge  =  0.0;
3636   Short_t  maxpos4    = -1;
3637   Float_t  maxcharge4 =  0.0;
3638   Short_t  maxpos5    = -1;
3639   Float_t  maxcharge5 =  0.0;
3640
3641   //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3642   //if (tchi2s[bestiter]>25.) sigma2=1000.;  // dont'accept
3643
3644   Double_t vD          = AliTRDcalibDB::Instance()->GetVdrift(0,0,0);
3645   Double_t exB         = AliTRDcalibDB::Instance()->GetOmegaTau(vD);
3646   Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
3647   if (mpads > 3.5) {
3648     expectederr +=  (mpads - 3.5) * 0.04;
3649   }
3650   if (changes[bestiter] > 1) {
3651     expectederr +=  changes[bestiter] * 0.01;
3652   }
3653   expectederr += (0.03 * (tany-exB)*(tany-exB)) * 15.0;
3654   //if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
3655   //expectederr+=10000;
3656
3657   for (Int_t it = 0; it < t1-t0; it++) {
3658
3659     if (!cl[best[bestiter][it]][it]) continue;
3660     // Set cluster error
3661     cl[best[bestiter][it]][it]->SetSigmaY2(expectederr);
3662     if (!cl[best[bestiter][it]][it]->IsUsed()) {
3663       cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY());
3664       //cl[best[bestiter][it]][it]->Use();
3665     }
3666
3667     // Time bins with maximal charge
3668     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge ) {
3669       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3670       maxpos    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3671     }
3672     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3673       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
3674         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3675         maxpos4    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3676       }
3677     }
3678     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3679       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
3680         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3681         maxpos5    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3682       }
3683     }
3684
3685     // Time bins with maximal charge
3686     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge ) {
3687       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3688       maxpos    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3689     }
3690     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge4) {
3691       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 4) {
3692         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3693         maxpos4    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3694       }
3695     }
3696     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ()) > maxcharge5) {
3697       if (cl[best[bestiter][it]][it]->GetLocalTimeBin() >= 5) {
3698         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3699         maxpos5    = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3700       }
3701     }
3702
3703     clusters[it+t0] = indexes[best[bestiter][it]][it];
3704
3705   }
3706
3707   //
3708   // Set tracklet parameters
3709   //
3710   Double_t trackleterr2 = smoffset[bestiter] + 0.01*0.01;
3711   if (mpads > 3.5) {
3712     trackleterr2 += (mpads - 3.5) * 0.04;
3713   }
3714   trackleterr2 += changes[bestiter] * 0.01;
3715   trackleterr2 *= TMath::Max(14.0 - nfound,1.0);
3716   trackleterr2 += 0.2 * (tany-exB)*(tany-exB);
3717
3718   tracklet.Set(xmean
3719               ,track2.GetY()+moffset[bestiter]
3720               ,meanz[bestiter]
3721               ,track2.GetAlpha()
3722               ,trackleterr2);
3723   tracklet.SetTilt(h01);
3724   tracklet.SetP0(mean[bestiter]);
3725   tracklet.SetP1(angle[bestiter]);
3726   tracklet.SetN(nfound);
3727   tracklet.SetNCross(changes[bestiter]);
3728   tracklet.SetPlane(plane);
3729   tracklet.SetSigma2(expectederr);
3730   tracklet.SetChi2(tchi2s[bestiter]);
3731   tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
3732   track->fTracklets[plane]  = tracklet;
3733   track->fNWrong           += nbad[0];
3734
3735   //
3736   // Debuging part
3737   //
3738   TClonesArray array0("AliTRDcluster");
3739   TClonesArray array1("AliTRDcluster");
3740   array0.ExpandCreateFast(t1 - t0 + 1);
3741   array1.ExpandCreateFast(t1 - t0 + 1);
3742   TTreeSRedirector &cstream = *fDebugStreamer;
3743   AliTRDcluster dummy;
3744   Double_t dy0[100];
3745   Double_t dyb[100];
3746
3747   for (Int_t it = 0; it < t1-t0; it++) {
3748
3749     dy0[it] = dy[0][it];
3750     dyb[it] = dy[best[bestiter][it]][it];
3751     if(cl[0][it]) {
3752       new(array0[it]) AliTRDcluster(*cl[0][it]);
3753     }
3754     else{
3755       new(array0[it]) AliTRDcluster(dummy);
3756     }
3757     if(cl[best[bestiter][it]][it]) {
3758       new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
3759     }
3760     else{
3761       new(array1[it]) AliTRDcluster(dummy);
3762     }
3763
3764   }
3765
3766   TGraph graph0(t1-t0,x,dy0);
3767   TGraph graph1(t1-t0,x,dyb);
3768   TGraph graphy(t1-t0,x,yt);
3769   TGraph graphz(t1-t0,x,zt);
3770
3771   if (AliTRDReconstructor::StreamLevel() > 0) {
3772     cstream << "tracklet"
3773             << "track.="      << track              // Track parameters
3774             << "tany="        << tany               // Tangent of the local track angle
3775             << "xmean="       << xmean              // xmean - reference x of tracklet
3776             << "tilt="        << h01                // Tilt angle
3777             << "nall="        << nall               // Number of foundable clusters
3778             << "nfound="      << nfound             // Number of found clusters
3779             << "clfound="     << clfound            // Total number of found clusters in road
3780             << "mpads="       << mpads              // Mean number of pads per cluster
3781             << "plane="       << plane              // Plane number
3782             << "detector="    << detector           // Detector number
3783             << "road="        << road               // The width of the used road
3784             << "graph0.="     << &graph0            // x - y = dy for closest cluster
3785             << "graph1.="     << &graph1            // x - y = dy for second closest cluster
3786             << "graphy.="     << &graphy            // y position of the track
3787             << "graphz.="     << &graphz            // z position of the track
3788             //<< "fCl.="        << &array0            // Closest cluster
3789             //<< "fCl2.="       << &array1            // Second closest cluster
3790             << "maxpos="      << maxpos             // Maximal charge postion
3791             << "maxcharge="   << maxcharge          // Maximal charge
3792             << "maxpos4="     << maxpos4            // Maximal charge postion - after bin 4
3793             << "maxcharge4="  << maxcharge4         // Maximal charge         - after bin 4
3794             << "maxpos5="     << maxpos5            // Maximal charge postion - after bin 5
3795             << "maxcharge5="  << maxcharge5         // Maximal charge         - after bin 5
3796             << "bestiter="    << bestiter           // Best iteration number
3797             << "tracklet.="   << &tracklet          // Corrspond to the best iteration
3798             << "tchi20="      << tchi2s[0]          // Chi2 of cluster in the 0 iteration
3799             << "tchi2b="      << tchi2s[bestiter]   // Chi2 of cluster in the best  iteration
3800             << "sigmas0="     << sigmas[0]          // Residuals sigma
3801             << "sigmasb="     << sigmas[bestiter]   // Residulas sigma
3802             << "ngood0="      << ngood[0]           // Number of good clusters in 0 iteration
3803             << "nbad0="       << nbad[0]            // Number of bad clusters in 0 iteration
3804             << "ngoodb="      << ngood[bestiter]    // Number of bad clusters in best iteration
3805             << "nbadb="       << nbad[bestiter]     // Number of good clusters in best iteration
3806             << "changes0="    << changes[0]         // Changes of pardrows in iteration number 0
3807             << "changesb="    << changes[bestiter]  // Changes of pardrows in best iteration
3808             << "moffset0="    << moffset[0]         // Offset fixing angle in iter=0
3809             << "smoffset0="   << smoffset[0]        // Sigma of offset fixing angle in iter=0
3810             << "moffsetb="    << moffset[bestiter]  // Offset fixing angle in iter=best
3811             << "smoffsetb="   << smoffset[bestiter] // Sigma of offset fixing angle in iter=best
3812             << "mean0="       << mean[0]            // Mean dy in iter=0;
3813             << "smean0="      << smean[0]           // Sigma of mean dy in iter=0
3814             << "meanb="       << mean[bestiter]     // Mean dy in iter=best
3815             << "smeanb="      << smean[bestiter]    // Sigma of mean dy in iter=best
3816             << "angle0="      << angle[0]           // Angle deviation in the iteration number 0
3817             << "sangle0="     << sangle[0]          // Sigma of angular deviation in iter=0
3818             << "angleb="      << angle[bestiter]    // Angle deviation in the best iteration
3819             << "sangleb="     << sangle[bestiter]   // Sigma of angle deviation in the best iteration
3820             << "expectederr=" << expectederr        // Expected error of cluster position
3821             << "\n";
3822   }
3823
3824   return nfound;
3825
3826 }
3827
3828 //_____________________________________________________________________________
3829 Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
3830                         , Int_t *outlist, Bool_t down)
3831 {
3832   //
3833   //  Sort eleements according occurancy
3834   //  The size of output array has is 2*n
3835   //
3836
3837   Int_t *sindexS = new Int_t[n];     // Temp array for sorting
3838   Int_t *sindexF = new Int_t[2*n];
3839   for (Int_t i = 0; i < n; i++) {
3840     sindexF[i] = 0;
3841   }
3842
3843   TMath::Sort(n,inlist,sindexS,down);
3844   Int_t last     = inlist[sindexS[0]];
3845   Int_t val      = last;
3846   sindexF[0]     = 1;
3847   sindexF[0+n]   = last;
3848   Int_t countPos = 0;
3849
3850   // Find frequency
3851   for (Int_t i = 1; i < n; i++) {
3852     val = inlist[sindexS[i]];
3853     if (last == val) {
3854       sindexF[countPos]++;
3855     }
3856     else {
3857       countPos++;
3858       sindexF[countPos+n] = val;
3859       sindexF[countPos]++;
3860       last                = val;
3861     }
3862   }
3863   if (last == val) countPos++;
3864
3865   // Sort according frequency
3866   TMath::Sort(countPos,sindexF,sindexS,kTRUE);
3867   for (Int_t i = 0; i < countPos; i++) {
3868     outlist[2*i  ] = sindexF[sindexS[i]+n];
3869     outlist[2*i+1] = sindexF[sindexS[i]];
3870   }
3871
3872   delete [] sindexS;
3873   delete [] sindexF;
3874
3875   return countPos;
3876
3877 }
3878
3879 //_____________________________________________________________________________
3880 AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
3881 {
3882   //
3883   // Register a seed
3884   //
3885
3886   Double_t alpha = AliTRDgeometry::GetAlpha();
3887   Double_t shift = AliTRDgeometry::GetAlpha() / 2.0;
3888
3889   Double_t c[15];
3890   c[ 0] = 0.2;
3891   c[ 1] = 0.0; c[ 2] = 2.0;
3892   c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
3893   c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0;  c[ 9] = 0.1;
3894   c[10] = 0.0; c[11] = 0.0; c[12] = 0.0;  c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
3895
3896   Int_t          index = 0;
3897   AliTRDcluster *cl    = 0;
3898
3899   for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
3900     if (seeds[ilayer].IsOK()) {
3901       for (Int_t itime = 22; itime > 0; itime--) {
3902         if (seeds[ilayer].fIndexes[itime] > 0) {
3903           index = seeds[ilayer].fIndexes[itime];
3904           cl    = seeds[ilayer].fClusters[itime];
3905           break;
3906         }
3907       }
3908     }
3909     if (index > 0) break;
3910   }
3911   if (cl == 0) return 0;
3912
3913   AliTRDtrack *track = new AliTRDtrack(cl
3914                                       ,index
3915                                       ,&params[1]
3916                                       ,c
3917                                       ,params[0]
3918                                       ,params[6]*alpha+shift);
3919   track->PropagateTo(params[0]-5.0);
3920   track->ResetCovariance(1);
3921
3922   Int_t rc = FollowBackProlongation(*track);
3923   if (rc < 30) {
3924     delete track;
3925     track = 0;
3926   }
3927   else {
3928     track->CookdEdx();
3929     CookdEdxTimBin(*track);
3930     CookLabel(track,0.9);
3931   }
3932
3933   return track;
3934
3935 }