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