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