TRD/AliTRDtracker.cxx

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