TRD/AliTRDtracker.cxx

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