TRD/AliTRDtracker.cxx

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