TRD/AliTRDtrackerV1.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id$ */
  18
  19 ///////////////////////////////////////////////////////////////////////////////
  20 //                                                                           //
  21 //  Track finder                                                             //
  22 //                                                                           //
  23 //  Authors:                                                                 //
  24 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
  25 //    Markus Fasel <M.Fasel@gsi.de>                                          //
  26 //                                                                           //
  27 ///////////////////////////////////////////////////////////////////////////////
  28
  29 #include <Riostream.h>
  30 #include <stdio.h>
  31 #include <string.h>
  32
  33 #include <TBranch.h>
  34 #include <TFile.h>
  35 #include <TGraph.h>
  36 #include <TH1D.h>
  37 #include <TH2D.h>
  38 #include <TLinearFitter.h>
  39 #include <TROOT.h>
  40 #include <TTree.h>
  41 #include <TClonesArray.h>
  42 #include <TRandom.h>
  43 #include <TTreeStream.h>
  44
  45 #include "AliLog.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAlignObj.h"
  48 #include "AliRieman.h"
  49 #include "AliTrackPointArray.h"
  50
  51 #include "AliTRDtrackerV1.h"
  52 #include "AliTRDtrackingChamber.h"
  53 #include "AliTRDgeometry.h"
  54 #include "AliTRDpadPlane.h"
  55 #include "AliTRDgeometry.h"
  56 #include "AliTRDcluster.h"
  57 #include "AliTRDtrack.h"
  58 #include "AliTRDseed.h"
  59 #include "AliTRDcalibDB.h"
  60 #include "AliTRDCommonParam.h"
  61 #include "AliTRDReconstructor.h"
  62 #include "AliTRDCalibraFillHisto.h"
  63 #include "AliTRDchamberTimeBin.h"
  64 #include "AliTRDrecoParam.h"
  65 #include "AliTRDseedV1.h"
  66 #include "AliTRDtrackV1.h"
  67 #include "Cal/AliTRDCalDet.h"
  68
  69
  70 ClassImp(AliTRDtrackerV1)
  71
  72
  73 const  Float_t  AliTRDtrackerV1::fgkMinClustersInTrack =  0.5;  //
  74 const  Float_t  AliTRDtrackerV1::fgkLabelFraction      =  0.8;  //
  75 const  Double_t AliTRDtrackerV1::fgkMaxChi2            = 12.0;  //
  76 const  Double_t AliTRDtrackerV1::fgkMaxSnp             =  0.95; // Maximum local sine of the azimuthal angle
  77 const  Double_t AliTRDtrackerV1::fgkMaxStep            =  2.0;  // Maximal step size in propagation
  78 Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
  79                 0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
  80                 0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
  81                 0.0474, 0.0408, 0.0335, 0.0335, 0.0335
  82 };
  83 TTreeSRedirector *AliTRDtrackerV1::fgDebugStreamer = 0x0;
  84 AliRieman* AliTRDtrackerV1::fgRieman = 0x0;
  85 TLinearFitter* AliTRDtrackerV1::fgTiltedRieman = 0x0;
  86 TLinearFitter* AliTRDtrackerV1::fgTiltedRiemanConstrained = 0x0;
  87
  88 //____________________________________________________________________
  89 AliTRDtrackerV1::AliTRDtrackerV1()
  90   :AliTracker()
  91   ,fGeom(new AliTRDgeometry())
  92   ,fClusters(0x0)
  93   ,fTracklets(0x0)
  94   ,fTracks(0x0)
  95   ,fTimeBinsPerPlane(0)
  96   ,fSieveSeeding(0)
  97 {
  98   //
  99   // Default constructor.
 100   //
 101   if (!AliTRDcalibDB::Instance()) {
 102     AliFatal("Could not get calibration object");
 103   }
 104   fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
 105
 106         for (Int_t isector = 0; isector < AliTRDgeometry::kNsect; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector, fTimeBinsPerPlane);
 107
 108   if(AliTRDReconstructor::StreamLevel() > 1){
 109                 TDirectory *savedir = gDirectory;
 110                 fgDebugStreamer    = new TTreeSRedirector("TRD.TrackerDebug.root");
 111         savedir->cd();
 112         }
 113 }
 114
 115 //____________________________________________________________________
 116 AliTRDtrackerV1::~AliTRDtrackerV1()
 117 {
 118   //
 119   // Destructor
 120   //
 121
 122         if(fgDebugStreamer) delete fgDebugStreamer;
 123         if(fgRieman) delete fgRieman;
 124         if(fgTiltedRieman) delete fgTiltedRieman;
 125         if(fgTiltedRiemanConstrained) delete fgTiltedRiemanConstrained;
 126         if(fTracks) {fTracks->Delete(); delete fTracks;}
 127         if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
 128         if(fClusters) {fClusters->Delete(); delete fClusters;}
 129         if(fGeom) delete fGeom;
 130 }
 131
 132 //____________________________________________________________________
 133 Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
 134 {
 135   //
 136   // Steering stand alone tracking for full TRD detector
 137   //
 138   // Parameters :
 139   //   esd     : The ESD event. On output it contains
 140   //             the ESD tracks found in TRD.
 141   //
 142   // Output :
 143   //   Number of tracks found in the TRD detector.
 144   //
 145   // Detailed description
 146   // 1. Launch individual SM trackers.
 147   //    See AliTRDtrackerV1::Clusters2TracksSM() for details.
 148   //
 149
 150         if(!AliTRDReconstructor::RecoParam()){
 151                 AliError("Reconstruction configuration not initialized. Call first AliTRDReconstructor::SetRecoParam().");
 152                 return 0;
 153         }
 154
 155         //AliInfo("Start Track Finder ...");
 156         Int_t ntracks = 0;
 157         for(int ism=0; ism<AliTRDgeometry::kNsect; ism++){
 158 //      for(int ism=1; ism<2; ism++){
 159                         //AliInfo(Form("Processing supermodule %i ...", ism));
 160                         ntracks += Clusters2TracksSM(ism, esd);
 161         }
 162   AliInfo(Form("Number of found tracks : %d", ntracks));
 163         return ntracks;
 164 }
 165
 166
 167 //_____________________________________________________________________________
 168 Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t index, AliTrackPoint &p) const
 169 {
 170         //AliInfo(Form("Asking for tracklet %d", index));
 171
 172         if(index<0) return kFALSE;
 173         AliTRDseedV1 *tracklet = 0x0;
 174         if(!(tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index))) return kFALSE;
 175
 176         // get detector for this tracklet
 177         AliTRDcluster *cl = 0x0;
 178         Int_t ic = 0; do; while(!(cl = tracklet->GetClusters(ic++)));
 179         Int_t  idet     = cl->GetDetector();
 180
 181         Double_t local[3];
 182         local[0] = tracklet->GetX0();
 183         local[1] = tracklet->GetYfit(0);
 184         local[2] = tracklet->GetZfit(0);
 185         Double_t global[3];
 186         fGeom->RotateBack(idet, local, global);
 187         p.SetXYZ(global[0],global[1],global[2]);
 188
 189
 190         // setting volume id
 191         AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
 192         switch (fGeom->GetPlane(idet)) {
 193         case 0:
 194                 iLayer = AliGeomManager::kTRD1;
 195                 break;
 196         case 1:
 197                 iLayer = AliGeomManager::kTRD2;
 198                 break;
 199         case 2:
 200                 iLayer = AliGeomManager::kTRD3;
 201                 break;
 202         case 3:
 203                 iLayer = AliGeomManager::kTRD4;
 204                 break;
 205         case 4:
 206                 iLayer = AliGeomManager::kTRD5;
 207                 break;
 208         case 5:
 209                 iLayer = AliGeomManager::kTRD6;
 210                 break;
 211         };
 212         Int_t    modId = fGeom->GetSector(idet) * fGeom->Ncham() + fGeom->GetChamber(idet);
 213         UShort_t volid = AliGeomManager::LayerToVolUID(iLayer, modId);
 214         p.SetVolumeID(volid);
 215
 216         return kTRUE;
 217 }
 218
 219 //____________________________________________________________________
 220 TLinearFitter* AliTRDtrackerV1::GetTiltedRiemanFitter()
 221 {
 222         if(!fgTiltedRieman) fgTiltedRieman = new TLinearFitter(4, "hyp4");
 223         return fgTiltedRieman;
 224 }
 225
 226 //____________________________________________________________________
 227 TLinearFitter* AliTRDtrackerV1::GetTiltedRiemanFitterConstraint()
 228 {
 229         if(!fgTiltedRiemanConstrained) fgTiltedRiemanConstrained = new TLinearFitter(2, "hyp2");
 230         return fgTiltedRiemanConstrained;
 231 }
 232
 233 //____________________________________________________________________
 234 AliRieman* AliTRDtrackerV1::GetRiemanFitter()
 235 {
 236         if(!fgRieman) fgRieman = new AliRieman(AliTRDtrackingChamber::kNTimeBins * AliTRDgeometry::kNplan);
 237         return fgRieman;
 238 }
 239
 240 //_____________________________________________________________________________
 241 Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
 242 {
 243   //
 244   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
 245   // backpropagated by the TPC tracker. Each seed is first propagated
 246   // to the TRD, and then its prolongation is searched in the TRD.
 247   // If sufficiently long continuation of the track is found in the TRD
 248   // the track is updated, otherwise it's stored as originaly defined
 249   // by the TPC tracker.
 250   //
 251
 252         Int_t   found    = 0;     // number of tracks found
 253         Float_t foundMin = 20.0;
 254
 255         AliTRDseed::SetNTimeBins(fTimeBinsPerPlane);
 256         Int_t    nSeed   = event->GetNumberOfTracks();
 257         if(!nSeed){
 258                 // run stand alone tracking
 259                 if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
 260                 return 0;
 261         }
 262
 263         Float_t *quality = new Float_t[nSeed];
 264         Int_t   *index   = new Int_t[nSeed];
 265         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
 266                 AliESDtrack *seed = event->GetTrack(iSeed);
 267                 Double_t covariance[15];
 268                 seed->GetExternalCovariance(covariance);
 269                 quality[iSeed] = covariance[0] + covariance[2];
 270         }
 271         // Sort tracks according to covariance of local Y and Z
 272         TMath::Sort(nSeed,quality,index,kFALSE);
 273
 274         // Backpropagate all seeds
 275         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
 276
 277                 // Get the seeds in sorted sequence
 278                 AliESDtrack *seed = event->GetTrack(index[iSeed]);
 279
 280                 // Check the seed status
 281                 ULong_t status = seed->GetStatus();
 282                 if ((status & AliESDtrack::kTPCout) == 0) continue;
 283                 if ((status & AliESDtrack::kTRDout) != 0) continue;
 284
 285                 // Do the back prolongation
 286                 Int_t   lbl         = seed->GetLabel();
 287                 AliTRDtrackV1 *track  = new AliTRDtrackV1(*seed);
 288                 //track->Print();
 289                 track->SetSeedLabel(lbl);
 290                 seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); // Make backup
 291                 Float_t p4          = track->GetC();
 292                 Int_t   expectedClr = FollowBackProlongation(*track);
 293                 //AliInfo(Form("\nTRACK %d Clusters %d [%d] in chi2 %f", index[iSeed], expectedClr, track->GetNumberOfClusters(), track->GetChi2()));
 294                 //track->Print();
 295
 296                 //Double_t cov[15];
 297                 //seed->GetExternalCovariance(cov);
 298                 //AliInfo(Form("track %d cov[%f %f] 0", index[iSeed], cov[0], cov[2]));
 299
 300                 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
 301                                 (track->Pt() > 0.8)) {
 302                         //
 303                         // Make backup for back propagation
 304                         //
 305                         Int_t foundClr = track->GetNumberOfClusters();
 306                         if (foundClr >= foundMin) {
 307                                 //AliInfo(Form("Making backup track ncls [%d]...", foundClr));
 308                                 track->CookdEdx();
 309                                 track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
 310                                 CookLabel(track,1 - fgkLabelFraction);
 311                                 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
 312
 313
 314                 //seed->GetExternalCovariance(cov);
 315                 //AliInfo(Form("track %d cov[%f %f] 0 test", index[iSeed], cov[0], cov[2]));
 316
 317                                 // Sign only gold tracks
 318                                 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
 319                                         if ((seed->GetKinkIndex(0)      ==   0) &&
 320                                                         (track->Pt()                <  1.5)) UseClusters(track);
 321                                 }
 322                                 Bool_t isGold = kFALSE;
 323
 324                                 // Full gold track
 325                                 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
 326                                         if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 327
 328                                         isGold = kTRUE;
 329                                 }
 330                 //seed->GetExternalCovariance(cov);
 331                 //AliInfo(Form("track %d cov[%f %f] 00", index[iSeed], cov[0], cov[2]));
 332
 333                                 // Almost gold track
 334                                 if ((!isGold)  && (track->GetNCross() == 0) &&
 335                                                 (track->GetChi2() / track->GetNumberOfClusters()  < 7)) {
 336                                         //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 337                                         if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 338
 339                                         isGold = kTRUE;
 340                                 }
 341                 //seed->GetExternalCovariance(cov);
 342                 //AliInfo(Form("track %d cov[%f %f] 01", index[iSeed], cov[0], cov[2]));
 343
 344                                 if ((!isGold) && (track->GetBackupTrack())) {
 345                                         if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
 346                                                 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 347                                                 isGold = kTRUE;
 348                                         }
 349                                 }
 350                 //seed->GetExternalCovariance(cov);
 351                 //AliInfo(Form("track %d cov[%f %f] 02", index[iSeed], cov[0], cov[2]));
 352
 353                                 //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected())  > 0.4)) {
 354                                         //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 355                                 //}
 356                         }
 357                 }
 358                 /**/
 359
 360                 /**/
 361                 // Debug part of tracking
 362 /*              TTreeSRedirector &cstream = *fgDebugStreamer;
 363                 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.
 364                 if (AliTRDReconstructor::StreamLevel() > 0) {
 365                         if (track->GetBackupTrack()) {
 366                                 cstream << "Tracks"
 367                                 << "EventNrInFile="  << eventNrInFile
 368                                 << "ESD.="     << seed
 369                                 << "trd.="     << track
 370                                 << "trdback.=" << track->GetBackupTrack()
 371                                 << "\n";
 372                         }
 373                         else {
 374                                 cstream << "Tracks"
 375                                 << "EventNrInFile="  << eventNrInFile
 376                                 << "ESD.="     << seed
 377                                 << "trd.="     << track
 378                                 << "trdback.=" << track
 379                                 << "\n";
 380                         }
 381                 }*/
 382                 /**/
 383
 384                 //seed->GetExternalCovariance(cov);
 385                 //AliInfo(Form("track %d cov[%f %f] 1", index[iSeed], cov[0], cov[2]));
 386
 387                 // Propagation to the TOF (I.Belikov)
 388                 if (track->GetStop() == kFALSE) {
 389                         //AliInfo("Track not stopped in TRD ...");
 390                         Double_t xtof  = 371.0;
 391                         Double_t xTOF0 = 370.0;
 392
 393                         Double_t c2    = track->GetSnp() + track->GetC() * (xtof - track->GetX());
 394                         if (TMath::Abs(c2) >= 0.99) {
 395                                 delete track;
 396                                 continue;
 397                         }
 398
 399                         PropagateToX(*track,xTOF0,fgkMaxStep);
 400
 401                         // Energy losses taken to the account - check one more time
 402                         c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
 403                         if (TMath::Abs(c2) >= 0.99) {
 404                                 delete track;
 405                                 continue;
 406                         }
 407
 408                         //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
 409                         //      fHBackfit->Fill(7);
 410                         //delete track;
 411                         //      continue;
 412                         //}
 413
 414                         Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
 415                         Double_t y;
 416                         track->GetYAt(xtof,GetBz(),y);
 417                         if (y >  ymax) {
 418                                 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
 419                                         delete track;
 420                                         continue;
 421                                 }
 422                         }else if (y < -ymax) {
 423                                 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
 424                                         delete track;
 425                                         continue;
 426                                 }
 427                         }
 428
 429                         if (track->PropagateTo(xtof)) {
 430                                 //AliInfo("set kTRDout");
 431                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
 432
 433                                 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 434                                         for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
 435                                                 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
 436                                         }
 437                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 438                                 }
 439                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
 440                                 if (track->GetNumberOfClusters() > foundMin) found++;
 441                         }
 442                 } else {
 443                         //AliInfo("Track stopped in TRD ...");
 444
 445                         if ((track->GetNumberOfClusters() >              15) &&
 446                                         (track->GetNumberOfClusters() > 0.5*expectedClr)) {
 447                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
 448
 449                                 //seed->SetStatus(AliESDtrack::kTRDStop);
 450                                 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 451                                         for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
 452                                                 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
 453                                         }
 454                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 455                                 }
 456                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
 457                                 found++;
 458                         }
 459                 }
 460
 461                 //if (((t->GetStatus()&AliESDtrack::kTRDout)!=0 )
 462
 463                 seed->SetTRDQuality(track->StatusForTOF());
 464                 seed->SetTRDBudget(track->GetBudget(0));
 465                 delete track;
 466         }
 467
 468
 469         AliInfo(Form("Number of seeds: %d", nSeed));
 470         AliInfo(Form("Number of back propagated TRD tracks: %d", found));
 471
 472         delete [] index;
 473         delete [] quality;
 474
 475   return 0;
 476 }
 477
 478
 479 //____________________________________________________________________
 480 Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event)
 481 {
 482   //
 483   // Refits tracks within the TRD. The ESD event is expected to contain seeds
 484   // at the outer part of the TRD.
 485   // The tracks are propagated to the innermost time bin
 486   // of the TRD and the ESD event is updated
 487   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
 488   //
 489
 490   Int_t   nseed    = 0; // contor for loaded seeds
 491   Int_t   found    = 0; // contor for updated TRD tracks
 492
 493   // Calibration monitor
 494   AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
 495   if (!calibra) AliInfo("Could not get Calibra instance\n");
 496
 497
 498   AliTRDtrackV1 track;
 499   for (Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++) {
 500     AliESDtrack *seed = event->GetTrack(itrack);
 501                 new(&track) AliTRDtrackV1(*seed);
 502
 503     if (track.GetX() < 270.0) {
 504       seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup);
 505       //AliInfo(Form("Remove for X = %7.3f [270.]\n", track.GetX()));
 506                         continue;
 507     }
 508
 509     ULong_t status = seed->GetStatus();
 510     if((status & AliESDtrack::kTRDout) == 0) continue;
 511     if((status & AliESDtrack::kTRDin)  != 0) continue;
 512     nseed++;
 513
 514     track.ResetCovariance(50.0);
 515
 516                 // do the propagation and processing
 517     Bool_t kUPDATE = kFALSE;
 518                 Double_t xTPC = 250.0;
 519     if(FollowProlongation(track)){
 520                         // computes PID for track
 521                         track.CookPID();
 522                         // update calibration references using this track
 523                         if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(&track);
 524
 525                         // Prolongate to TPC
 526                         if (PropagateToX(track, xTPC, fgkMaxStep)) { //  -with update
 527                                 seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit);
 528                                 track.UpdateESDtrack(seed);
 529                                 // Add TRD track to ESDfriendTrack
 530                                 if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
 531                                         AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
 532                                         calibTrack->SetOwner();
 533                                         seed->AddCalibObject(calibTrack);
 534                                 }
 535                                 found++;
 536                                 kUPDATE = kTRUE;
 537                         }
 538                 }
 539
 540                 // Prolongate to TPC without update
 541                 if(!kUPDATE) {
 542       AliTRDtrackV1 tt(*seed);
 543       if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDrefit);
 544     }
 545   }
 546   AliInfo(Form("Number of loaded seeds: %d",nseed));
 547   AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
 548
 549         return 0;
 550 }
 551
 552
 553 //____________________________________________________________________
 554 Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
 555 {
 556 // Extrapolates the TRD track in the TPC direction.
 557 //
 558 // Parameters
 559 //   t : the TRD track which has to be extrapolated
 560 //
 561 // Output
 562 //   number of clusters attached to the track
 563 //
 564 // Detailed description
 565 //
 566 // Starting from current radial position of track <t> this function
 567 // extrapolates the track through the 6 TRD layers. The following steps
 568 // are being performed for each plane:
 569 // 1. prepare track:
 570 //   a. get plane limits in the local x direction
 571 //   b. check crossing sectors
 572 //   c. check track inclination
 573 // 2. search tracklet in the tracker list (see GetTracklet() for details)
 574 // 3. evaluate material budget using the geo manager
 575 // 4. propagate and update track using the tracklet information.
 576 //
 577 // Debug level 2
 578 //
 579
 580         //AliInfo("");
 581         Int_t    nClustersExpected = 0;
 582         Int_t lastplane = 5; //GetLastPlane(&t);
 583         for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
 584     Int_t   index   = 0;
 585     AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
 586                 if(!tracklet) continue;
 587                 if(!tracklet->IsOK()) AliWarning("tracklet not OK");
 588
 589                 t.SetTracklet(tracklet, iplane, index);
 590
 591                 Double_t x  = tracklet->GetX0();
 592     if (x < (t.GetX()-fgkMaxStep) && !PropagateToX(t, x+fgkMaxStep, fgkMaxStep)) break;
 593     if (!AdjustSector(&t)) break;
 594
 595     // Start global position
 596     Double_t xyz0[3];
 597     t.GetXYZ(xyz0);
 598
 599                 // End global position
 600     Double_t alpha = t.GetAlpha(), y, z;
 601     if (!t.GetProlongation(x,y,z)) break;
 602     Double_t xyz1[3];
 603     xyz1[0] =  x * TMath::Cos(alpha) - y * TMath::Sin(alpha);
 604     xyz1[1] =  x * TMath::Sin(alpha) + y * TMath::Cos(alpha);
 605     xyz1[2] =  z;
 606
 607     // Get material budget
 608     Double_t param[7];
 609     AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
 610     Double_t xrho= param[0]*param[4];
 611     Double_t xx0 = param[1]; // Get mean propagation parameters
 612
 613     // Propagate and update
 614                 t.PropagateTo(x, xx0, xrho);
 615           if (!AdjustSector(&t)) break;
 616
 617     Double_t maxChi2 = t.GetPredictedChi2(tracklet);
 618           if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){
 619                 nClustersExpected += tracklet->GetN();
 620         }
 621   }
 622
 623         if(AliTRDReconstructor::StreamLevel() > 1){
 624                 Int_t index;
 625                 for(int iplane=0; iplane<6; iplane++){
 626                         AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
 627                         if(!tracklet) continue;
 628                         t.SetTracklet(tracklet, iplane, index);
 629                 }
 630
 631                 TTreeSRedirector &cstreamer = *fgDebugStreamer;
 632                 cstreamer << "FollowProlongation"
 633                         << "ncl="      << nClustersExpected
 634                         << "track.="   << &t
 635                         << "\n";
 636         }
 637
 638   return nClustersExpected;
 639
 640 }
 641
 642 //_____________________________________________________________________________
 643 Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
 644 {
 645 // Extrapolates the TRD track in the TOF direction.
 646 //
 647 // Parameters
 648 //   t : the TRD track which has to be extrapolated
 649 //
 650 // Output
 651 //   number of clusters attached to the track
 652 //
 653 // Detailed description
 654 //
 655 // Starting from current radial position of track <t> this function
 656 // extrapolates the track through the 6 TRD layers. The following steps
 657 // are being performed for each plane:
 658 // 1. prepare track:
 659 //   a. get plane limits in the local x direction
 660 //   b. check crossing sectors
 661 //   c. check track inclination
 662 // 2. build tracklet (see AliTRDseed::AttachClusters() for details)
 663 // 3. evaluate material budget using the geo manager
 664 // 4. propagate and update track using the tracklet information.
 665 //
 666 // Debug level 2
 667 //
 668
 669         Int_t nClustersExpected = 0;
 670   Double_t clength = AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick();
 671   AliTRDtrackingChamber *chamber = 0x0;
 672
 673   // Loop through the TRD planes
 674   for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
 675                 // BUILD TRACKLET IF NOT ALREADY BUILT
 676                 Double_t x = 0., y, z, alpha;
 677     AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
 678                 if(!tracklet.IsOK()){
 679                 alpha = t.GetAlpha();
 680                 Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsect));
 681
 682                 if(!fTrSec[sector].GetNChambers()) continue;
 683
 684                 if((x = fTrSec[sector].GetX(iplane)) < 1.) continue;
 685
 686                         if (!t.GetProlongation(x, y, z)) break;
 687                         Int_t stack = fGeom->GetChamber(z, iplane);
 688                         Int_t nCandidates = stack >= 0 ? 1 : 2;
 689                         z -= stack >= 0 ? 0. : 4.;
 690
 691                         for(int icham=0; icham<nCandidates; icham++, z+=8){
 692                                 if((stack = fGeom->GetChamber(z, iplane)) < 0) continue;
 693
 694                                 if(!(chamber = fTrSec[sector].GetChamber(stack, iplane))) continue;
 695
 696                                 if(chamber->GetNClusters() < fTimeBinsPerPlane*AliTRDReconstructor::AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
 697
 698                                 x = chamber->GetX();
 699
 700                                 AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, stack);
 701                                 tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
 702                                 tracklet.SetPadLength(pp->GetLengthIPad());
 703                                 tracklet.SetPlane(iplane);
 704                                 tracklet.SetX0(x);
 705                                 tracklet.Init(&t);
 706                                 if(!tracklet.AttachClustersIter(chamber, 1000.)) continue;
 707                                 tracklet.Init(&t);
 708
 709                 if(tracklet.GetN() < fTimeBinsPerPlane * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
 710
 711                                 break;
 712                         }
 713                 }
 714     if(!tracklet.IsOK()){
 715                         if(x < 1.) continue; //temporary
 716                         if(!PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) break;
 717                         if(!AdjustSector(&t)) break;
 718                         if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
 719         continue;
 720     }
 721
 722                 // Propagate closer to the current chamber if neccessary
 723     x -= clength;
 724     if (x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) break;
 725     if (!AdjustSector(&t)) break;
 726     if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
 727
 728                 // load tracklet to the tracker and the track
 729                 Int_t index = SetTracklet(&tracklet);
 730                 t.SetTracklet(&tracklet, iplane, index);
 731
 732
 733                 // Calculate the mean material budget along the path inside the chamber
 734     //Calculate global entry and exit positions of the track in chamber (only track prolongation)
 735     Double_t xyz0[3]; // entry point
 736                 t.GetXYZ(xyz0);
 737                 alpha = t.GetAlpha();
 738                 x = tracklet.GetX0();
 739                 if (!t.GetProlongation(x, y, z)) break;
 740                 Double_t xyz1[3]; // exit point
 741                 xyz1[0] =  x * TMath::Cos(alpha) - y * TMath::Sin(alpha);
 742     xyz1[1] = +x * TMath::Sin(alpha) + y * TMath::Cos(alpha);
 743     xyz1[2] =  z;
 744     Double_t param[7];
 745                 AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
 746     // The mean propagation parameters
 747     Double_t xrho = param[0]*param[4]; // density*length
 748     Double_t xx0  = param[1]; // radiation length
 749
 750                 // Propagate and update track
 751                 t.PropagateTo(x, xx0, xrho);
 752           if (!AdjustSector(&t)) break;
 753                 Double_t maxChi2 = t.GetPredictedChi2(&tracklet);
 754                 if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){
 755                         nClustersExpected += tracklet.GetN();
 756                 }
 757                 // Reset material budget if 2 consecutive gold
 758                 if(iplane>0 && tracklet.GetN() + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
 759
 760                 // Make backup of the track until is gold
 761                 // TO DO update quality check of the track.
 762                 // consider comparison with fTimeBinsRange
 763                 Float_t ratio0 = tracklet.GetN() / Float_t(fTimeBinsPerPlane);
 764                 //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
 765     //printf("tracklet.GetChi2() %f     [< 18.0]\n", tracklet.GetChi2());
 766                 //printf("ratio0    %f              [>   0.8]\n", ratio0);
 767                 //printf("ratio1     %f             [>   0.6]\n", ratio1);
 768                 //printf("ratio0+ratio1 %f          [>   1.5]\n", ratio0+ratio1);
 769                 //printf("t.GetNCross()  %d         [==    0]\n", t.GetNCross());
 770                 //printf("TMath::Abs(t.GetSnp()) %f [<  0.85]\n", TMath::Abs(t.GetSnp()));
 771                 //printf("t.GetNumberOfClusters() %d [>    20]\n", t.GetNumberOfClusters());
 772
 773                 if (//(tracklet.GetChi2()      <  18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update
 774         (ratio0                  >   0.8) &&
 775         //(ratio1                  >   0.6) &&
 776         //(ratio0+ratio1           >   1.5) &&
 777         (t.GetNCross()           ==    0) &&
 778         (TMath::Abs(t.GetSnp())  <  0.85) &&
 779         (t.GetNumberOfClusters() >    20)) t.MakeBackupTrack();
 780
 781         } // end planes loop
 782
 783         if(AliTRDReconstructor::StreamLevel() > 1){
 784                 TTreeSRedirector &cstreamer = *fgDebugStreamer;
 785                 cstreamer << "FollowBackProlongation"
 786                         << "ncl="      << nClustersExpected
 787                         << "track.="   << &t
 788                         << "\n";
 789         }
 790
 791         return nClustersExpected;
 792 }
 793
 794 //_________________________________________________________________________
 795 Float_t AliTRDtrackerV1::FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_t *planes){
 796 //
 797 // Fits a Riemann-circle to the given points without tilting pad correction.
 798 // The fit is performed using an instance of the class AliRieman (equations
 799 // and transformations see documentation of this class)
 800 // Afterwards all the tracklets are Updated
 801 //
 802 // Parameters: - Array of tracklets (AliTRDseedV1)
 803 //             - Storage for the chi2 values (beginning with direction z)
 804 //             - Seeding configuration
 805 // Output:     - The curvature
 806 //
 807   AliRieman *fitter = AliTRDtrackerV1::GetRiemanFitter();
 808         fitter->Reset();
 809   Int_t allplanes[] = {0, 1, 2, 3, 4, 5};
 810   Int_t *ppl = &allplanes[0];
 811         Int_t maxLayers = 6;
 812   if(planes){
 813     maxLayers = 4;
 814     ppl = planes;
 815   }
 816   for(Int_t il = 0; il < maxLayers; il++){
 817                 if(!tracklets[ppl[il]].IsOK()) continue;
 818     fitter->AddPoint(tracklets[ppl[il]].GetX0(), tracklets[ppl[il]].GetYfitR(0), tracklets[ppl[il]].GetZProb(),1,10);
 819   }
 820   fitter->Update();
 821   // Set the reference position of the fit and calculate the chi2 values
 822   memset(chi2, 0, sizeof(Double_t) * 2);
 823         for(Int_t il = 0; il < maxLayers; il++){
 824                 // Reference positions
 825                 tracklets[ppl[il]].Init(fitter);
 826
 827                 // chi2
 828                 if((!tracklets[ppl[il]].IsOK()) && (!planes)) continue;
 829                 chi2[0] += tracklets[ppl[il]].GetChi2Z();
 830                 chi2[1] += tracklets[ppl[il]].GetChi2Y();
 831         }
 832         return fitter->GetC();
 833 }
 834
 835 //_________________________________________________________________________
 836 void AliTRDtrackerV1::FitRieman(AliTRDcluster **seedcl, Double_t chi2[2])
 837 {
 838 //
 839 // Performs a Riemann helix fit using the seedclusters as spacepoints
 840 // Afterwards the chi2 values are calculated and the seeds are updated
 841 //
 842 // Parameters: - The four seedclusters
 843 //             - The tracklet array (AliTRDseedV1)
 844 //             - The seeding configuration
 845 //             - Chi2 array
 846 //
 847 // debug level 2
 848 //
 849         AliRieman *fitter = AliTRDtrackerV1::GetRiemanFitter();
 850         fitter->Reset();
 851         for(Int_t i = 0; i < 4; i++)
 852                 fitter->AddPoint(seedcl[i]->GetX(), seedcl[i]->GetY(), seedcl[i]->GetZ(), 1, 10);
 853         fitter->Update();
 854
 855
 856         // Update the seed and calculated the chi2 value
 857         chi2[0] = 0; chi2[1] = 0;
 858         for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
 859                 // chi2
 860                 chi2[0] += (seedcl[ipl]->GetZ() - fitter->GetZat(seedcl[ipl]->GetX())) * (seedcl[ipl]->GetZ() - fitter->GetZat(seedcl[ipl]->GetX()));
 861                 chi2[1] += (seedcl[ipl]->GetY() - fitter->GetYat(seedcl[ipl]->GetX())) * (seedcl[ipl]->GetY() - fitter->GetYat(seedcl[ipl]->GetX()));
 862         }
 863 }
 864
 865
 866 //_________________________________________________________________________
 867 Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Double_t zVertex)
 868 {
 869 //
 870 // Fits a helix to the clusters. Pad tilting is considered. As constraint it is
 871 // assumed that the vertex position is set to 0.
 872 // This method is very usefull for high-pt particles
 873 // Basis for the fit: (x - x0)^2 + (y - y0)^2 - R^2 = 0
 874 //      x0, y0: Center of the circle
 875 // Measured y-position: ymeas = y - tan(phiT)(zc - zt)
 876 //      zc: center of the pad row
 877 // Equation which has to be fitted (after transformation):
 878 // a + b * u + e * v + 2*(ymeas + tan(phiT)(z - zVertex))*t = 0
 879 // Transformation:
 880 // t = 1/(x^2 + y^2)
 881 // u = 2 * x * t
 882 // v = 2 * x * tan(phiT) * t
 883 // Parameters in the equation:
 884 //    a = -1/y0, b = x0/y0, e = dz/dx
 885 //
 886 // The Curvature is calculated by the following equation:
 887 //               - curv = a/Sqrt(b^2 + 1) = 1/R
 888 // Parameters:   - the 6 tracklets
 889 //               - the Vertex constraint
 890 // Output:       - the Chi2 value of the track
 891 //
 892 // debug level 5
 893 //
 894
 895         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 896         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 897
 898         TLinearFitter *fitter = GetTiltedRiemanFitterConstraint();
 899         fitter->StoreData(kTRUE);
 900         fitter->ClearPoints();
 901
 902         Float_t x, y, z, w, t, error, tilt;
 903         Double_t uvt[2];
 904         Int_t nPoints = 0;
 905   for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
 906                 if(!tracklets[ipl].IsOK()) continue;
 907                 for(Int_t itb = 0; itb < nTimeBins; itb++){
 908                         if(!tracklets[ipl].IsUsable(itb)) continue;
 909                         x = tracklets[ipl].GetX(itb) + tracklets[ipl].GetX0();
 910                         y = tracklets[ipl].GetY(itb);
 911                         z = tracklets[ipl].GetZ(itb);
 912                         tilt = tracklets[ipl].GetTilt();
 913                         // Transformation
 914                         t = 1/(x * x + y * y);
 915                         uvt[0] = 2 * x* t;
 916                         uvt[1] =  2.0 * tilt * x * t;
 917                         w = 2.0 * (y + tilt * (z - zVertex)) * t;
 918                         error = 2 * 0.2 * t;
 919                         fitter->AddPoint(uvt, w, error);
 920                         nPoints++;
 921                 }
 922         }
 923         fitter->Eval();
 924
 925         // Calculate curvature
 926         Double_t a = fitter->GetParameter(0);
 927         Double_t b = fitter->GetParameter(0);
 928         Double_t curvature = a/TMath::Sqrt(b*b + 1);
 929
 930         Float_t chi2track = fitter->GetChisquare()/Double_t(nPoints);
 931         for(Int_t ip = 0; ip < AliTRDtrackerV1::kNPlanes; ip++)
 932                 tracklets[ip].SetCC(curvature);
 933
 934         if(AliTRDReconstructor::StreamLevel() >= 5){
 935                 //Linear Model on z-direction
 936           Double_t xref = (tracklets[2].GetX0() + tracklets[3].GetX0())/2;              // Relative to the middle of the stack
 937                 Double_t slope = fitter->GetParameter(2);
 938                 Double_t zref = slope * xref;
 939                 Float_t chi2Z = CalculateChi2Z(tracklets, zref, slope);
 940                 TTreeSRedirector &treeStreamer = *fgDebugStreamer;
 941                 treeStreamer << "FitTiltedRiemanConstraint"
 942                         << "Curvature=" << curvature
 943                         << "Chi2Track=" << chi2track
 944                         << "Chi2Z="                     << chi2Z
 945                         << "zref="                      << zref
 946                         << "\n";
 947         }
 948         return chi2track;
 949 }
 950
 951 //_________________________________________________________________________
 952 Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigError)
 953 {
 954 //
 955 // Performs a Riemann fit taking tilting pad correction into account
 956 // The equation of a Riemann circle, where the y position is substituted by the
 957 // measured y-position taking pad tilting into account, has to be transformed
 958 // into a 4-dimensional hyperplane equation
 959 // Riemann circle: (x-x0)^2 + (y-y0)^2 -R^2 = 0
 960 // Measured y-Position: ymeas = y - tan(phiT)(zc - zt)
 961 //          zc: center of the pad row
 962 //          zt: z-position of the track
 963 // The z-position of the track is assumed to be linear dependent on the x-position
 964 // Transformed equation: a + b * u + c * t + d * v  + e * w - 2 * (ymeas + tan(phiT) * zc) * t = 0
 965 // Transformation:       u = 2 * x * t
 966 //                       v = 2 * tan(phiT) * t
 967 //                       w = 2 * tan(phiT) * (x - xref) * t
 968 //                       t = 1 / (x^2 + ymeas^2)
 969 // Parameters:           a = -1/y0
 970 //                       b = x0/y0
 971 //                       c = (R^2 -x0^2 - y0^2)/y0
 972 //                       d = offset
 973 //                       e = dz/dx
 974 // If the offset respectively the slope in z-position is impossible, the parameters are fixed using
 975 // results from the simple riemann fit. Afterwards the fit is redone.
 976 // The curvature is calculated according to the formula:
 977 //                       curv = a/(1 + b^2 + c*a) = 1/R
 978 //
 979 // Paramters:   - Array of tracklets (connected to the track candidate)
 980 //              - Flag selecting the error definition
 981 // Output:      - Chi2 value of the track
 982 //
 983 // debug level 5
 984 //
 985
 986         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 987         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 988
 989         TLinearFitter *fitter = GetTiltedRiemanFitter();
 990   fitter->StoreData(kTRUE);
 991         fitter->ClearPoints();
 992
 993         // Calculate the reference position:
 994         Int_t nDistances = 0;
 995         Float_t meanDistance = 0.;
 996         Int_t startIndex = 5;
 997         for(Int_t il =5; il > 0; il--){
 998                 if(tracklets[il].IsOK() && tracklets[il -1].IsOK()){
 999                         meanDistance += tracklets[il].GetX0() - tracklets[il -1].GetX0();
1000                         nDistances++;
1001                 }
1002                 if(tracklets[il].IsOK()) startIndex = il;
1003         }
1004         meanDistance /= nDistances;
1005         if(tracklets[0].IsOK()) startIndex = 0;
1006         Double_t xref = tracklets[startIndex].GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
1007
1008         Float_t x, y, z, t, tilt, xdelta, rhs, error;
1009         Float_t dzMean = 0;     Int_t dzcounter = 0;    // A reference z and a reference slope is used if the fitresults in z-direction are not acceptable
1010         Double_t uvt[4];
1011         Int_t nPoints = 0;
1012         for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
1013                 if(!tracklets[ipl].IsOK()) continue;
1014                 dzMean += tracklets[ipl].GetZfitR(1);
1015                 dzcounter++;
1016                 for(Int_t itb = 0; itb < nTimeBins; itb++){
1017                         if (!tracklets[ipl].IsUsable(itb)) continue;
1018                         x = tracklets[ipl].GetX(itb) + tracklets[ipl].GetX0();
1019                         y = tracklets[ipl].GetY(itb);
1020                         z = tracklets[ipl].GetZ(itb);
1021                         tilt = tracklets[ipl].GetTilt();
1022                         xdelta = x - xref;
1023                         // Transformation
1024                         t = 1/(x*x + y*y);
1025                         uvt[0] = 2.0 * x * t;
1026                         uvt[1] = t;
1027                         uvt[2] = 2.0 * tilt * t;
1028                         uvt[3] = 2.0 * tilt * xdelta * t;
1029                         rhs = 2.0 * (y + tilt*z) * t;
1030                         // error definition changes for the different calls
1031                         error = 2.0 * t;
1032                         error *= sigError ? tracklets[ipl].GetSigmaY() : 0.2;
1033                         fitter->AddPoint(uvt, rhs, error);
1034                         nPoints++;
1035                 }
1036         }
1037
1038         fitter->Eval();
1039
1040         Double_t offset = fitter->GetParameter(3);
1041         Double_t slope  = fitter->GetParameter(4);
1042
1043         // Linear fitter  - not possible to make boundaries
1044         // Do not accept non possible z and dzdx combinations
1045         Bool_t acceptablez = kTRUE;
1046         Double_t zref = 0.0;
1047         for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNplan; iLayer++) {
1048                 if(!tracklets[iLayer].IsOK()) continue;
1049                 zref = offset + slope * (tracklets[iLayer].GetX0() - xref);
1050                 if (TMath::Abs(tracklets[iLayer].GetZProb() - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0)
1051                         acceptablez = kFALSE;
1052         }
1053         if (!acceptablez) {
1054                 dzMean /= dzcounter;
1055                 Double_t zmf  = tracklets[startIndex].GetZfitR(0) + dzMean * (xref - tracklets[startIndex].GetX0()); // Z-Position of the track at the middle of a stack assuming a linear dependence on x (approximation)
1056                 fgTiltedRieman->FixParameter(3, zmf);
1057                 fgTiltedRieman->FixParameter(4, dzMean);
1058                 fitter->Eval();
1059                 fitter->ReleaseParameter(3);
1060                 fitter->ReleaseParameter(4);
1061                 offset = fitter->GetParameter(3);
1062                 slope = fitter->GetParameter(4);
1063         }
1064
1065         // Calculate Curvarture
1066         Double_t a     =  fitter->GetParameter(0);
1067         Double_t b     =  fitter->GetParameter(1);
1068         Double_t c     =  fitter->GetParameter(2);
1069         Double_t curvature =  1.0 + b*b - c*a;
1070         Double_t dca  =  0.0;                                                                                                                   // Distance to closest approach
1071         if (curvature > 0.0) {
1072                 dca = -c / (TMath::Sqrt(1.0 + b*b - c*a) + TMath::Sqrt(1.0 + b*b));
1073                 curvature  =  a / TMath::Sqrt(curvature);
1074         }
1075
1076         Double_t chi2track = fitter->GetChisquare()/Double_t(nPoints);
1077
1078         // Update the tracklets
1079         Double_t dy, dz;
1080         for(Int_t iLayer = 0; iLayer < AliTRDtrackerV1::kNPlanes; iLayer++) {
1081
1082                 x  = tracklets[iLayer].GetX0();
1083                 y  = 0;
1084                 z  = 0;
1085                 dy = 0;
1086                 dz = 0;
1087
1088                 // y:     R^2 = (x - x0)^2 + (y - y0)^2
1089                 //     =>   y = y0 +/- Sqrt(R^2 - (x - x0)^2)
1090                 //          R = Sqrt() = 1/Curvature
1091                 //     =>   y = y0 +/- Sqrt(1/Curvature^2 - (x - x0)^2)
1092                 Double_t res = (x * a + b);                                                             // = (x - x0)/y0
1093                 res *= res;
1094                 res  = 1.0 - c * a + b * b - res;                                       // = (R^2 - (x - x0)^2)/y0^2
1095                 if (res >= 0) {
1096                         res = TMath::Sqrt(res);
1097                         y    = (1.0 - res) / a;
1098                 }
1099
1100                 // dy:      R^2 = (x - x0)^2 + (y - y0)^2
1101                 //     =>     y = +/- Sqrt(R^2 - (x - x0)^2) + y0
1102                 //     => dy/dx = (x - x0)/Sqrt(R^2 - (x - x0)^2)
1103                 // Curvature: cr = 1/R = a/Sqrt(1 + b^2 - c*a)
1104                 //     => dy/dx =  (x - x0)/(1/(cr^2) - (x - x0)^2)
1105                 Double_t x0 = -b / a;
1106                 if (-c * a + b * b + 1 > 0) {
1107                         if (1.0/(curvature * curvature) - (x - x0) * (x - x0) > 0.0) {
1108                                 Double_t yderiv = (x - x0) / TMath::Sqrt(1.0/(curvature * curvature) - (x - x0) * (x - x0));
1109                                 if (a < 0) yderiv *= -1.0;
1110                                 dy = yderiv;
1111                         }
1112                 }
1113                 z  = offset + slope * (x - xref);
1114                 dz = slope;
1115                 tracklets[iLayer].SetYref(0, y);
1116                 tracklets[iLayer].SetYref(1, dy);
1117                 tracklets[iLayer].SetZref(0, z);
1118                 tracklets[iLayer].SetZref(1, dz);
1119                 tracklets[iLayer].SetC(curvature);
1120                 tracklets[iLayer].SetChi2(chi2track);
1121   }
1122
1123
1124         if(AliTRDReconstructor::StreamLevel() >= 5){
1125                 Double_t chi2Z = CalculateChi2Z(tracklets, offset, slope);
1126                 TTreeSRedirector &treeStreamer = *fgDebugStreamer;
1127                 treeStreamer << "FitTiltedRieman"
1128                         << "error="         << sigError
1129                         << "Curvature="     << curvature
1130                         << "Chi2track="     << chi2track
1131                         << "Chi2Z="         << chi2Z
1132                         << "D="             << c
1133                         << "DCA="           << dca
1134                         << "Offset="        << offset
1135                         << "Slope="         << slope
1136                         << "\n";
1137         }
1138
1139         return chi2track;
1140 }
1141
1142 //_________________________________________________________________________
1143 Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope)
1144 {
1145 //
1146 // Calculates the chi2-value of the track in z-Direction including tilting pad correction.
1147 // A linear dependence on the x-value serves as a model.
1148 // The parameters are related to the tilted Riemann fit.
1149 // Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate
1150 //             - the offset for the reference x
1151 //             - the slope
1152 // Output:     - The Chi2 value of the track in z-Direction
1153 //
1154         Double_t xref = .5 * (tracklets[2].GetX0() + tracklets[3].GetX0());
1155         Float_t chi2Z = 0, nLayers = 0;
1156         for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNplan; iLayer++) {
1157                 if(!tracklets[iLayer].IsOK()) continue;
1158                 Double_t z = offset + slope * (tracklets[iLayer].GetX0() - xref);
1159                 chi2Z += TMath::Abs(tracklets[iLayer].GetMeanz() - z);
1160                 nLayers++;
1161         }
1162         chi2Z /= TMath::Max((nLayers - 3.0),1.0);
1163         return chi2Z;
1164 }
1165
1166
1167
1168 //_____________________________________________________________________________
1169 Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep)
1170 {
1171   //
1172   // Starting from current X-position of track <t> this function
1173   // extrapolates the track up to radial position <xToGo>.
1174   // Returns 1 if track reaches the plane, and 0 otherwise
1175   //
1176
1177   const Double_t kEpsilon = 0.00001;
1178
1179   // Current track X-position
1180   Double_t xpos = t.GetX();
1181
1182   // Direction: inward or outward
1183   Double_t dir  = (xpos < xToGo) ? 1.0 : -1.0;
1184
1185   while (((xToGo - xpos) * dir) > kEpsilon) {
1186
1187     Double_t xyz0[3];
1188     Double_t xyz1[3];
1189     Double_t param[7];
1190     Double_t x;
1191     Double_t y;
1192     Double_t z;
1193
1194     // The next step size
1195     Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
1196
1197     // Get the global position of the starting point
1198     t.GetXYZ(xyz0);
1199
1200     // X-position after next step
1201     x = xpos + step;
1202
1203     // Get local Y and Z at the X-position of the next step
1204     if (!t.GetProlongation(x,y,z)) {
1205       return 0; // No prolongation possible
1206     }
1207
1208     // The global position of the end point of this prolongation step
1209     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
1210     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
1211     xyz1[2] =  z;
1212
1213     // Calculate the mean material budget between start and
1214     // end point of this prolongation step
1215     AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
1216
1217     // Propagate the track to the X-position after the next step
1218     if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
1219       return 0;
1220     }
1221
1222     // Rotate the track if necessary
1223     AdjustSector(&t);
1224
1225     // New track X-position
1226     xpos = t.GetX();
1227
1228   }
1229
1230   return 1;
1231
1232 }
1233
1234
1235 //_____________________________________________________________________________
1236 Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) const
1237 {
1238   //
1239   // Reads AliTRDclusters from the file.
1240   // The names of the cluster tree and branches
1241   // should match the ones used in AliTRDclusterizer::WriteClusters()
1242   //
1243
1244   Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster)));
1245   TObjArray *clusterArray = new TObjArray(nsize+1000);
1246
1247   TBranch *branch = clusterTree->GetBranch("TRDcluster");
1248   if (!branch) {
1249     AliError("Can't get the branch !");
1250     return 1;
1251   }
1252   branch->SetAddress(&clusterArray);
1253
1254   if(!fClusters){
1255         array = new TClonesArray("AliTRDcluster", nsize);
1256         array->SetOwner(kTRUE);
1257   }
1258
1259   // Loop through all entries in the tree
1260   Int_t nEntries   = (Int_t) clusterTree->GetEntries();
1261   Int_t nbytes     = 0;
1262   Int_t ncl        = 0;
1263   AliTRDcluster *c = 0x0;
1264   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1265     // Import the tree
1266     nbytes += clusterTree->GetEvent(iEntry);
1267
1268     // Get the number of points in the detector
1269     Int_t nCluster = clusterArray->GetEntriesFast();
1270     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
1271       if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
1272       new((*fClusters)[ncl++]) AliTRDcluster(*c);
1273       clusterArray->RemoveAt(iCluster);
1274     }
1275
1276   }
1277   delete clusterArray;
1278
1279   return 0;
1280 }
1281
1282 //_____________________________________________________________________________
1283 Int_t AliTRDtrackerV1::LoadClusters(TTree *cTree)
1284 {
1285   //
1286   // Fills clusters into TRD tracking_sectors
1287   // Note that the numbering scheme for the TRD tracking_sectors
1288   // differs from that of TRD sectors
1289   //
1290
1291
1292   if (ReadClusters(fClusters, cTree)) {
1293     AliError("Problem with reading the clusters !");
1294     return 1;
1295   }
1296   Int_t ncl  = fClusters->GetEntriesFast(), nin = 0;
1297   Int_t icl = ncl;
1298   while (icl--) {
1299     AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(icl);
1300                 if(c->IsInChamber()) nin++;
1301     Int_t detector       = c->GetDetector();
1302     Int_t sector         = fGeom->GetSector(detector);
1303     Int_t stack          = fGeom->GetChamber(detector);
1304     Int_t plane          = fGeom->GetPlane(detector);
1305
1306                 fTrSec[sector].GetChamber(stack, plane, kTRUE)->InsertCluster(c, icl);
1307   }
1308   AliInfo(Form("Clusters %d in %6.2f %%", ncl, 100.*float(nin)/ncl));
1309
1310         for(int isector =0; isector<AliTRDgeometry::kNsect; isector++){
1311                 if(!fTrSec[isector].GetNChambers()) continue;
1312                 fTrSec[isector].Init();
1313   }
1314
1315   return 0;
1316 }
1317
1318
1319 //____________________________________________________________________
1320 void AliTRDtrackerV1::UnloadClusters()
1321 {
1322   //
1323   // Clears the arrays of clusters and tracks. Resets sectors and timebins
1324   //
1325
1326         if(fTracks) fTracks->Delete();
1327   if(fTracklets) fTracklets->Delete();
1328   if(fClusters) fClusters->Delete();
1329
1330   for (int i = 0; i < AliTRDgeometry::kNsect; i++) fTrSec[i].Clear();
1331
1332 }
1333
1334 //_____________________________________________________________________________
1335 Bool_t AliTRDtrackerV1::AdjustSector(AliTRDtrackV1 *track)
1336 {
1337   //
1338   // Rotates the track when necessary
1339   //
1340
1341   Double_t alpha = AliTRDgeometry::GetAlpha();
1342   Double_t y     = track->GetY();
1343   Double_t ymax  = track->GetX()*TMath::Tan(0.5*alpha);
1344
1345   if      (y >  ymax) {
1346     if (!track->Rotate( alpha)) {
1347       return kFALSE;
1348     }
1349   }
1350   else if (y < -ymax) {
1351     if (!track->Rotate(-alpha)) {
1352       return kFALSE;
1353     }
1354   }
1355
1356   return kTRUE;
1357
1358 }
1359
1360
1361 //____________________________________________________________________
1362 AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t &idx)
1363 {
1364 // Find tracklet for TRD track <track>
1365 // Parameters
1366 // - track
1367 // - sector
1368 // - plane
1369 // - index
1370 // Output
1371 // tracklet
1372 // index
1373 // Detailed description
1374 //
1375         idx = track->GetTrackletIndex(p);
1376         AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
1377
1378         return tracklet;
1379 }
1380
1381 //____________________________________________________________________
1382 Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet)
1383 {
1384 // Add this tracklet to the list of tracklets stored in the tracker
1385 //
1386 // Parameters
1387 //   - tracklet : pointer to the tracklet to be added to the list
1388 //
1389 // Output
1390 //   - the index of the new tracklet in the tracker tracklets list
1391 //
1392 // Detailed description
1393 // Build the tracklets list if it is not yet created (late initialization)
1394 // and adds the new tracklet to the list.
1395 //
1396         if(!fTracklets){
1397                 fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsect()*kMaxTracksStack);
1398                 fTracklets->SetOwner(kTRUE);
1399         }
1400         Int_t nentries = fTracklets->GetEntriesFast();
1401         new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
1402         return nentries;
1403 }
1404
1405 //____________________________________________________________________
1406 Int_t AliTRDtrackerV1::Clusters2TracksSM(Int_t sector, AliESDEvent *esd)
1407 {
1408   //
1409   // Steer tracking for one SM.
1410   //
1411   // Parameters :
1412   //   sector  : Array of (SM) propagation layers containing clusters
1413   //   esd     : The current ESD event. On output it contains the also
1414   //             the ESD (TRD) tracks found in this SM.
1415   //
1416   // Output :
1417   //   Number of tracks found in this TRD supermodule.
1418   //
1419   // Detailed description
1420   //
1421   // 1. Unpack AliTRDpropagationLayers objects for each stack.
1422   // 2. Launch stack tracking.
1423   //    See AliTRDtrackerV1::Clusters2TracksStack() for details.
1424   // 3. Pack results in the ESD event.
1425   //
1426
1427         // allocate space for esd tracks in this SM
1428         TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
1429         esdTrackList.SetOwner();
1430
1431         Int_t nTracks   = 0;
1432         Int_t nChambers = 0;
1433         AliTRDtrackingChamber **stack = 0x0, *chamber = 0x0;
1434         for(int istack = 0; istack<AliTRDgeometry::kNcham; istack++){
1435                 if(!(stack = fTrSec[sector].GetStack(istack))) continue;
1436                 nChambers = 0;
1437                 for(int iplane=0; iplane<AliTRDgeometry::kNplan; iplane++){
1438                         if(!(chamber = stack[iplane])) continue;
1439                         if(chamber->GetNClusters() < fTimeBinsPerPlane * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
1440                         nChambers++;
1441                         //AliInfo(Form("sector %d stack %d plane %d clusters %d", sector, istack, iplane, chamber->GetNClusters()));
1442                 }
1443                 if(nChambers < 4) continue;
1444                 //AliInfo(Form("Doing stack %d", istack));
1445                 nTracks += Clusters2TracksStack(stack, &esdTrackList);
1446         }
1447         //AliInfo(Form("Found %d tracks in SM %d [%d]\n", nTracks, sector, esd->GetNumberOfTracks()));
1448
1449         for(int itrack=0; itrack<nTracks; itrack++)
1450           esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
1451
1452         return nTracks;
1453 }
1454
1455 //____________________________________________________________________
1456 Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList)
1457 {
1458   //
1459   // Make tracks in one TRD stack.
1460   //
1461   // Parameters :
1462   //   layer  : Array of stack propagation layers containing clusters
1463   //   esdTrackList  : Array of ESD tracks found by the stand alone tracker.
1464   //                   On exit the tracks found in this stack are appended.
1465   //
1466   // Output :
1467   //   Number of tracks found in this stack.
1468   //
1469   // Detailed description
1470   //
1471   // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
1472   // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations.
1473   //    See AliTRDtrackerV1::MakeSeeds() for more details.
1474   // 3. Arrange track candidates in decreasing order of their quality
1475   // 4. Classify tracks in 5 categories according to:
1476   //    a) number of layers crossed
1477   //    b) track quality
1478   // 5. Sign clusters by tracks in decreasing order of track quality
1479   // 6. Build AliTRDtrack out of seeding tracklets
1480   // 7. Cook MC label
1481   // 8. Build ESD track and register it to the output list
1482   //
1483
1484         AliTRDtrackingChamber *chamber = 0x0;
1485         AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
1486         Int_t pars[4]; // MakeSeeds parameters
1487
1488         //Double_t alpha = AliTRDgeometry::GetAlpha();
1489         //Double_t shift = .5 * alpha;
1490         Int_t configs[kNConfigs];
1491
1492         // Build initial seeding configurations
1493         Double_t quality = BuildSeedingConfigs(stack, configs);
1494         if(AliTRDReconstructor::StreamLevel() > 1){
1495         AliInfo(Form("Plane config %d %d %d Quality %f"
1496     , configs[0], configs[1], configs[2], quality));
1497         }
1498
1499         // Initialize contors
1500         Int_t ntracks,      // number of TRD track candidates
1501               ntracks1,     // number of registered TRD tracks/iter
1502               ntracks2 = 0; // number of all registered TRD tracks in stack
1503         fSieveSeeding = 0;
1504         do{
1505                 // Loop over seeding configurations
1506                 ntracks = 0; ntracks1 = 0;
1507                 for (Int_t iconf = 0; iconf<3; iconf++) {
1508                         pars[0] = configs[iconf];
1509                         pars[1] = ntracks;
1510                         ntracks = MakeSeeds(stack, &sseed[6*ntracks], pars);
1511                         if(ntracks == kMaxTracksStack) break;
1512                 }
1513                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
1514
1515                 if(!ntracks) break;
1516
1517                 // Sort the seeds according to their quality
1518                 Int_t sort[kMaxTracksStack];
1519                 TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
1520
1521                 // Initialize number of tracks so far and logic switches
1522                 Int_t ntracks0 = esdTrackList->GetEntriesFast();
1523                 Bool_t signedTrack[kMaxTracksStack];
1524                 Bool_t fakeTrack[kMaxTracksStack];
1525                 for (Int_t i=0; i<ntracks; i++){
1526                         signedTrack[i] = kFALSE;
1527                         fakeTrack[i] = kFALSE;
1528                 }
1529                 //AliInfo("Selecting track candidates ...");
1530
1531                 // Sieve clusters in decreasing order of track quality
1532                 Double_t trackParams[7];
1533 //              AliTRDseedV1 *lseed = 0x0;
1534                 Int_t jSieve = 0, candidates;
1535                 do{
1536                         //AliInfo(Form("\t\tITER = %i ", jSieve));
1537
1538                         // Check track candidates
1539                         candidates = 0;
1540                         for (Int_t itrack = 0; itrack < ntracks; itrack++) {
1541                                 Int_t trackIndex = sort[itrack];
1542                                 if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
1543
1544
1545                                 // Calculate track parameters from tracklets seeds
1546                                 Int_t labelsall[1000];
1547                                 Int_t nlabelsall = 0;
1548                                 Int_t naccepted  = 0;
1549                                 Int_t ncl        = 0;
1550                                 Int_t nused      = 0;
1551                                 Int_t nlayers    = 0;
1552                                 Int_t findable   = 0;
1553                                 for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
1554                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
1555                                         if(!sseed[jseed].IsOK()) continue;
1556                                         if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15) findable++;
1557
1558                                         sseed[jseed].UpdateUsed();
1559                                         ncl   += sseed[jseed].GetN2();
1560                                         nused += sseed[jseed].GetNUsed();
1561                                         nlayers++;
1562
1563                                         // Cooking label
1564                                         for (Int_t itime = 0; itime < fTimeBinsPerPlane; itime++) {
1565                                                 if(!sseed[jseed].IsUsable(itime)) continue;
1566                                                 naccepted++;
1567                                                 Int_t tindex = 0, ilab = 0;
1568                                                 while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
1569                                                         labelsall[nlabelsall++] = tindex;
1570                                                         ilab++;
1571                                                 }
1572                                         }
1573                                 }
1574                                 // Filter duplicated tracks
1575                                 if (nused > 30){
1576                                         //printf("Skip %d nused %d\n", trackIndex, nused);
1577                                         fakeTrack[trackIndex] = kTRUE;
1578                                         continue;
1579                                 }
1580                                 if (Float_t(nused)/ncl >= .25){
1581                                         //printf("Skip %d nused/ncl >= .25\n", trackIndex);
1582                                         fakeTrack[trackIndex] = kTRUE;
1583                                         continue;
1584                                 }
1585
1586                                 // Classify tracks
1587                                 Bool_t skip = kFALSE;
1588                                 switch(jSieve){
1589                                 case 0:
1590                                         if(nlayers < 6) {skip = kTRUE; break;}
1591                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1592                                         break;
1593
1594                                 case 1:
1595                                         if(nlayers < findable){skip = kTRUE; break;}
1596                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
1597                                         break;
1598
1599                                 case 2:
1600                                         if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
1601                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
1602                                         break;
1603
1604                                 case 3:
1605                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1606                                         break;
1607
1608                                 case 4:
1609                                         if (nlayers == 3){skip = kTRUE; break;}
1610                                         //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
1611                                         break;
1612                                 }
1613                                 if(skip){
1614                                         candidates++;
1615                                         //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
1616                                         continue;
1617                                 }
1618                                 signedTrack[trackIndex] = kTRUE;
1619
1620
1621                                 // Build track label - what happens if measured data ???
1622                                 Int_t labels[1000];
1623                                 Int_t outlab[1000];
1624                                 Int_t nlab = 0;
1625                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1626                                         Int_t jseed = kNPlanes*trackIndex+iLayer;
1627                                         if(!sseed[jseed].IsOK()) continue;
1628                                         for(int ilab=0; ilab<2; ilab++){
1629                                                 if(sseed[jseed].GetLabels(ilab) < 0) continue;
1630                                                 labels[nlab] = sseed[jseed].GetLabels(ilab);
1631                                                 nlab++;
1632                                         }
1633                                 }
1634                                 Freq(nlab,labels,outlab,kFALSE);
1635                                 Int_t   label     = outlab[0];
1636                                 Int_t   frequency = outlab[1];
1637                                 Freq(nlabelsall,labelsall,outlab,kFALSE);
1638                                 Int_t   label1    = outlab[0];
1639                                 Int_t   label2    = outlab[2];
1640                                 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
1641
1642
1643                                 // Sign clusters
1644                                 AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
1645                                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1646                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
1647                                         if(!sseed[jseed].IsOK()) continue;
1648                                         if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
1649                                         sseed[jseed].UseClusters();
1650                                         if(!cl){
1651                                                 Int_t ic = 0;
1652                                                 while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
1653                                                 clusterIndex =  sseed[jseed].GetIndexes(ic);
1654                                         }
1655                                 }
1656                                 if(!cl) continue;
1657
1658
1659                                 // Build track parameters
1660                                 AliTRDseedV1 *lseed =&sseed[trackIndex*6];
1661                                 Int_t idx = 0;
1662                                 while(idx<3 && !lseed->IsOK()) {
1663                                         idx++;
1664                                         lseed++;
1665                                 }
1666                                 Double_t cR = lseed->GetC();
1667                                 trackParams[1] = lseed->GetYref(0);
1668                                 trackParams[2] = lseed->GetZref(0);
1669                                 trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
1670                                 trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
1671                                 trackParams[5] = cR;
1672                                 trackParams[0] = lseed->GetX0();
1673                                 Int_t ich = 0; while(!(chamber = stack[ich])) ich++;
1674                                 trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift;      // Supermodule*/
1675
1676                                 if(AliTRDReconstructor::StreamLevel() > 1){
1677                                         AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]));
1678
1679                                         Int_t nclusters = 0;
1680                                         AliTRDseedV1 *dseed[6];
1681                                         for(int is=0; is<6; is++){
1682                                                 dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is]);
1683                                                 dseed[is]->SetOwner();
1684                                                 nclusters += sseed[is].GetN2();
1685                                         }
1686                                         //Int_t eventNrInFile = esd->GetEventNumberInFile();
1687                                         //AliInfo(Form("Number of clusters %d.", nclusters));
1688                                         TTreeSRedirector &cstreamer = *fgDebugStreamer;
1689                                         cstreamer << "Clusters2TracksStack"
1690                                                 << "Iter="      << fSieveSeeding
1691                                                 << "Like="      << fTrackQuality[trackIndex]
1692                                                 << "S0.="       << dseed[0]
1693                                                 << "S1.="       << dseed[1]
1694                                                 << "S2.="       << dseed[2]
1695                                                 << "S3.="       << dseed[3]
1696                                                 << "S4.="       << dseed[4]
1697                                                 << "S5.="       << dseed[5]
1698                                                 << "p0=" << trackParams[0]
1699                                                 << "p1=" << trackParams[1]
1700                                                 << "p2=" << trackParams[2]
1701                                                 << "p3=" << trackParams[3]
1702                                                 << "p4=" << trackParams[4]
1703                                                 << "p5=" << trackParams[5]
1704                                                 << "p6=" << trackParams[6]
1705                                                 << "Label="     << label
1706                                                 << "Label1="    << label1
1707                                                 << "Label2="    << label2
1708                                                 << "FakeRatio=" << fakeratio
1709                                                 << "Freq="      << frequency
1710                                                 << "Ncl="       << ncl
1711                                                 << "NLayers="   << nlayers
1712                                                 << "Findable="  << findable
1713                                                 << "NUsed="     << nused
1714                                                 << "\n";
1715                                 }
1716
1717                                 AliTRDtrackV1 *track = MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
1718                                 if(!track){
1719                                         //AliWarning("Fail to build a TRD Track.");
1720                                         continue;
1721                                 }
1722                                 //AliInfo("End of MakeTrack()");
1723                                 AliESDtrack esdTrack;
1724                                 esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
1725                                 esdTrack.SetLabel(track->GetLabel());
1726                                 new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
1727                                 ntracks1++;
1728                         }
1729
1730                         jSieve++;
1731                 } while(jSieve<5 && candidates); // end track candidates sieve
1732                 if(!ntracks1) break;
1733
1734                 // increment counters
1735                 ntracks2 += ntracks1;
1736                 fSieveSeeding++;
1737
1738                 // Rebuild plane configurations and indices taking only unused clusters into account
1739                 quality = BuildSeedingConfigs(stack, configs);
1740                 if(quality < 1.E-7) break; //AliTRDReconstructor::RecoParam()->GetPlaneQualityThreshold()) break;
1741
1742                 for(Int_t ip = 0; ip < kNPlanes; ip++){
1743                         if(!(chamber = stack[ip])) continue;
1744                         chamber->Build(fGeom);//Indices(fSieveSeeding);
1745                 }
1746
1747                 if(AliTRDReconstructor::StreamLevel() > 1){
1748                         AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
1749                 }
1750         } while(fSieveSeeding<10); // end stack clusters sieve
1751
1752
1753
1754         //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
1755
1756         return ntracks2;
1757 }
1758
1759 //___________________________________________________________________
1760 Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int_t *configs)
1761 {
1762   //
1763   // Assign probabilities to chambers according to their
1764   // capability of producing seeds.
1765   //
1766   // Parameters :
1767   //
1768   //   layers : Array of stack propagation layers for all 6 chambers in one stack
1769   //   configs : On exit array of configuration indexes (see GetSeedingConfig()
1770   // for details) in the decreasing order of their seeding probabilities.
1771   //
1772   // Output :
1773   //
1774   //  Return top configuration quality
1775   //
1776   // Detailed description:
1777   //
1778   // To each chamber seeding configuration (see GetSeedingConfig() for
1779   // the list of all configurations) one defines 2 quality factors:
1780   //  - an apriori topological quality (see GetSeedingConfig() for details) and
1781   //  - a data quality based on the uniformity of the distribution of
1782   //    clusters over the x range (time bins population). See CookChamberQA() for details.
1783   // The overall chamber quality is given by the product of this 2 contributions.
1784   //
1785
1786         Double_t chamberQ[kNPlanes];
1787         AliTRDtrackingChamber *chamber = 0x0;
1788         for(int iplane=0; iplane<kNPlanes; iplane++){
1789                 if(!(chamber = stack[iplane])) continue;
1790                 chamberQ[iplane] = (chamber = stack[iplane]) ?  chamber->GetQuality(fTimeBinsPerPlane) : 0.;
1791         }
1792
1793         Double_t tconfig[kNConfigs];
1794         Int_t planes[4];
1795         for(int iconf=0; iconf<kNConfigs; iconf++){
1796                 GetSeedingConfig(iconf, planes);
1797                 tconfig[iconf] = fgTopologicQA[iconf];
1798                 for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQ[planes[iplane]];
1799         }
1800
1801         TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
1802 //      AliInfo(Form("q[%d] = %f", configs[0], tconfig[configs[0]]));
1803 //      AliInfo(Form("q[%d] = %f", configs[1], tconfig[configs[1]]));
1804 //      AliInfo(Form("q[%d] = %f", configs[2], tconfig[configs[2]]));
1805
1806         return tconfig[configs[0]];
1807 }
1808
1809 //____________________________________________________________________
1810 Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar)
1811 {
1812   //
1813   // Make tracklet seeds in the TRD stack.
1814   //
1815   // Parameters :
1816   //   layers : Array of stack propagation layers containing clusters
1817   //   sseed  : Array of empty tracklet seeds. On exit they are filled.
1818   //   ipar   : Control parameters:
1819   //       ipar[0] -> seeding chambers configuration
1820   //       ipar[1] -> stack index
1821   //       ipar[2] -> number of track candidates found so far
1822   //
1823   // Output :
1824   //   Number of tracks candidates found.
1825   //
1826   // Detailed description
1827   //
1828   // The following steps are performed:
1829   // 1. Select seeding layers from seeding chambers
1830   // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
1831   //   The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
1832   //   this order. The parameters controling the range of accepted clusters in
1833   //   layer 0, 1, and 2 are defined in AliTRDchamberTimeBin::BuildCond().
1834   // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
1835   // 4. Initialize seeding tracklets in the seeding chambers.
1836   // 5. Filter 0.
1837   //   Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
1838   //   Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
1839   // 6. Attach clusters to seeding tracklets and find linear approximation of
1840   //   the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
1841   //   clusters used by current seeds should not exceed ... (25).
1842   // 7. Filter 1.
1843   //   All 4 seeding tracklets should be correctly constructed (see
1844   //   AliTRDseedV1::AttachClustersIter())
1845   // 8. Helix fit of the seeding tracklets
1846   // 9. Filter 2.
1847   //   Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
1848   // 10. Extrapolation of the helix fit to the other 2 chambers:
1849   //    a) Initialization of extrapolation tracklet with fit parameters
1850   //    b) Helix fit of tracklets
1851   //    c) Attach clusters and linear interpolation to extrapolated tracklets
1852   //    d) Helix fit of tracklets
1853   // 11. Improve seeding tracklets quality by reassigning clusters.
1854   //      See AliTRDtrackerV1::ImproveSeedQuality() for details.
1855   // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
1856   // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
1857   // 14. Cooking labels for tracklets. Should be done only for MC
1858   // 15. Register seeds.
1859   //
1860
1861         AliTRDtrackingChamber *chamber = 0x0;
1862         AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
1863         AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
1864         Int_t ncl, mcl; // working variable for looping over clusters
1865         Int_t index[AliTRDchamberTimeBin::kMaxClustersLayer], jndex[AliTRDchamberTimeBin::kMaxClustersLayer];
1866         // chi2 storage
1867         // chi2[0] = tracklet chi2 on the Z direction
1868         // chi2[1] = tracklet chi2 on the R direction
1869         Double_t chi2[4];
1870
1871
1872         // this should be data member of AliTRDtrack
1873         Double_t seedQuality[kMaxTracksStack];
1874
1875         // unpack control parameters
1876         Int_t config  = ipar[0];
1877         Int_t ntracks = ipar[1];
1878         Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);
1879
1880         // Init chambers geometry
1881         Int_t ic = 0; while(!(chamber = stack[ic])) ic++;
1882         Int_t istack = fGeom->GetChamber(chamber->GetDetector());
1883         Double_t hL[kNPlanes];       // Tilting angle
1884         Float_t padlength[kNPlanes]; // pad lenghts
1885         AliTRDpadPlane *pp = 0x0;
1886         for(int iplane=0; iplane<kNPlanes; iplane++){
1887                 pp                = fGeom->GetPadPlane(iplane, istack);
1888                 hL[iplane]        = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
1889                 padlength[iplane] = pp->GetLengthIPad();
1890         }
1891
1892         if(AliTRDReconstructor::StreamLevel() > 1){
1893                 AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
1894         }
1895
1896         Int_t nlayers = 0;
1897         AliTRDchamberTimeBin *layer[] = {0x0, 0x0, 0x0, 0x0};
1898         for(int isl=0; isl<kNSeedPlanes; isl++){
1899                 if(!(chamber = stack[planes[isl]])) continue;
1900                 if(!(layer[isl] = chamber->GetSeedingLayer(fGeom))) continue;
1901                 nlayers++;
1902                 //AliInfo(Form("seeding plane %d clusters %d", planes[isl], Int_t(*layer[isl])));
1903         }
1904         if(nlayers < 4) return 0;
1905
1906
1907         // Start finding seeds
1908         Double_t cond0[4], cond1[4], cond2[4];
1909         Int_t icl = 0;
1910         while((c[3] = (*layer[3])[icl++])){
1911                 if(!c[3]) continue;
1912                 layer[0]->BuildCond(c[3], cond0, 0);
1913                 layer[0]->GetClusters(cond0, index, ncl);
1914                 //printf("Found c[3] candidates 0 %d\n", ncl);
1915                 Int_t jcl = 0;
1916                 while(jcl<ncl) {
1917                         c[0] = (*layer[0])[index[jcl++]];
1918                         if(!c[0]) continue;
1919                         Double_t dx    = c[3]->GetX() - c[0]->GetX();
1920                         Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
1921                         Double_t phi   = (c[3]->GetY() - c[0]->GetY())/dx;
1922                         layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
1923                         layer[1]->GetClusters(cond1, jndex, mcl);
1924                         //printf("Found c[0] candidates 1 %d\n", mcl);
1925
1926                         Int_t kcl = 0;
1927                         while(kcl<mcl) {
1928                                 c[1] = (*layer[1])[jndex[kcl++]];
1929                                 if(!c[1]) continue;
1930                                 layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
1931                                 c[2] = layer[2]->GetNearestCluster(cond2);
1932                                 //printf("Found c[1] candidate 2 %p\n", c[2]);
1933                                 if(!c[2]) continue;
1934
1935 //                              AliInfo("Seeding clusters found. Building seeds ...");
1936 //                              for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %6.3f, y = %6.3f, z = %6.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
1937
1938                                 for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
1939
1940                                 FitRieman(c, chi2);
1941
1942                                 AliTRDseedV1 *tseed = 0x0;
1943                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
1944                                         Int_t jLayer = planes[iLayer];
1945                                         tseed = &cseed[jLayer];
1946                                         tseed->SetPlane(jLayer);
1947                                         tseed->SetTilt(hL[jLayer]);
1948                                         tseed->SetPadLength(padlength[jLayer]);
1949                                         tseed->SetX0(stack[jLayer]->GetX());
1950                                         tseed->Init(GetRiemanFitter());
1951                                 }
1952
1953                                 Bool_t isFake = kFALSE;
1954                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1955                                         if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1956                                         if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1957                                         if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1958
1959                                         Float_t yref[4];
1960                                         for(int il=0; il<4; il++) yref[il] = cseed[planes[il]].GetYref(0);
1961                                         Int_t ll = c[3]->GetLabel(0);
1962                                         TTreeSRedirector &cs0 = *fgDebugStreamer;
1963                                                         cs0 << "MakeSeeds0"
1964                                                         <<"isFake=" << isFake
1965                                                         <<"label=" << ll
1966                                                         <<"chi2z=" << chi2[0]
1967                                                         <<"chi2y=" << chi2[1]
1968                                                         <<"yref0=" << yref[0]
1969                                                         <<"yref1=" << yref[1]
1970                                                         <<"yref2=" << yref[2]
1971                                                         <<"yref3=" << yref[3]
1972                                                         <<"c0.="   << c[0]
1973                                                         <<"c1.="   << c[1]
1974                                                         <<"c2.="   << c[2]
1975                                                         <<"c3.="   << c[3]
1976                                                         <<"\n";
1977                                 }
1978
1979                                 if(chi2[0] > AliTRDReconstructor::RecoParam()->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
1980                                         //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
1981                                         continue;
1982                                 }
1983                                 if(chi2[1] > AliTRDReconstructor::RecoParam()->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
1984                                         //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
1985                                         continue;
1986                                 }
1987                                 //AliInfo("Passed chi2 filter.");
1988
1989                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1990                                         Float_t minmax[2] = { -100.0,  100.0 };
1991                                         for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1992                                                 Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
1993                                                 if (max < minmax[1]) minmax[1] = max;
1994                                                 Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
1995                                                 if (min > minmax[0]) minmax[0] = min;
1996                                         }
1997                                         Double_t xpos[4];
1998                                         for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
1999                                         TTreeSRedirector &cstreamer = *fgDebugStreamer;
2000                                                         cstreamer << "MakeSeeds1"
2001                                                 << "isFake=" << isFake
2002                                                 << "config="   << config
2003                                                 << "Cl0.="   << c[0]
2004                                                 << "Cl1.="   << c[1]
2005                                                 << "Cl2.="   << c[2]
2006                                                 << "Cl3.="   << c[3]
2007                                                 << "X0="     << xpos[0] //layer[sLayer]->GetX()
2008                                                 << "X1="     << xpos[1] //layer[sLayer + 1]->GetX()
2009                                                 << "X2="     << xpos[2] //layer[sLayer + 2]->GetX()
2010                                                 << "X3="     << xpos[3] //layer[sLayer + 3]->GetX()
2011                                                 << "Y2exp="  << cond2[0]
2012                                                 << "Z2exp="  << cond2[1]
2013                                                 << "Chi2R="  << chi2[0]
2014                                                 << "Chi2Z="  << chi2[1]
2015                                                 << "Seed0.=" << &cseed[planes[0]]
2016                                                 << "Seed1.=" << &cseed[planes[1]]
2017                                                 << "Seed2.=" << &cseed[planes[2]]
2018                                                 << "Seed3.=" << &cseed[planes[3]]
2019                                                 << "Zmin="   << minmax[0]
2020                                                 << "Zmax="   << minmax[1]
2021                                                 << "\n" ;
2022                                 }
2023
2024                                 // try attaching clusters to tracklets
2025                                 Int_t nUsedCl = 0;
2026                                 Int_t nlayers = 0;
2027                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
2028                                         Int_t jLayer = planes[iLayer];
2029                                         if(!cseed[jLayer].AttachClustersIter(stack[jLayer], 5., kFALSE, c[iLayer])) continue;
2030                                         nUsedCl += cseed[jLayer].GetNUsed();
2031                                         if(nUsedCl > 25) break;
2032                                         nlayers++;
2033                                 }
2034                                 if(nlayers < kNSeedPlanes){
2035                                         //AliInfo(Form("Failed updating all seeds %d [%d].", nlayers, kNSeedPlanes));
2036                                         continue;
2037                                 }
2038                                 // fit tracklets and cook likelihood
2039                                 FitRieman(&cseed[0], chi2, &planes[0]);
2040                                 Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
2041                                 if (TMath::Log(1.E-9 + like) < AliTRDReconstructor::RecoParam()->GetTrackLikelihood()){
2042                                         //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
2043                                         continue;
2044                                 }
2045                                 //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
2046
2047
2048                                 // book preliminary results
2049                                 seedQuality[ntracks] = like;
2050                                 fSeedLayer[ntracks]  = config;/*sLayer;*/
2051
2052                                 // attach clusters to the extrapolation seeds
2053                                 Int_t lextrap[2];
2054                                 GetExtrapolationConfig(config, lextrap);
2055                                 Int_t nusedf   = 0; // debug value
2056                                 for(int iLayer=0; iLayer<2; iLayer++){
2057                                         Int_t jLayer = lextrap[iLayer];
2058                                         if(!(chamber = stack[jLayer])) continue;
2059
2060                                         // prepare extrapolated seed
2061                                         cseed[jLayer].Reset();
2062                                         cseed[jLayer].SetPlane(jLayer);
2063                                         cseed[jLayer].SetTilt(hL[jLayer]);
2064                                         cseed[jLayer].SetX0(chamber->GetX());
2065                                         cseed[jLayer].SetPadLength(padlength[jLayer]);
2066
2067                                         // fit extrapolated seed
2068                                         FitTiltedRieman(cseed, kTRUE);
2069                                         if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
2070                                         if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
2071                                         AliTRDseedV1 tseed = cseed[jLayer];
2072                                         if(!tseed.AttachClustersIter(chamber, 1000.)) continue;
2073                                         cseed[jLayer] = tseed;
2074                                         nusedf += cseed[jLayer].GetNUsed(); // debug value
2075                                 }
2076                                 FitTiltedRieman(cseed, kTRUE);
2077                                 //AliInfo("Extrapolation done.");
2078
2079                                 if(ImproveSeedQuality(stack, cseed) < 4) continue;
2080                                 //AliInfo("Improve seed quality done.");
2081
2082                                 // fit full track and cook likelihoods
2083                                 Double_t curv = FitRieman(&cseed[0], chi2);
2084                                 Double_t chi2ZF = chi2[0] / TMath::Max((nlayers - 3.), 1.);
2085                                 Double_t chi2RF = chi2[1] / TMath::Max((nlayers - 3.), 1.);
2086
2087                                 // do the final track fitting (Once with vertex constraint and once without vertex constraint)
2088                                 Double_t chi2Vals[3];
2089                                 chi2Vals[0] = FitTiltedRieman(&cseed[0], kFALSE);
2090                                 chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ());
2091                                 chi2Vals[2] = chi2ZF;
2092                                 fTrackQuality[ntracks] = CalculateTrackLikelihood(&cseed[0], &chi2Vals[0]);
2093                                 //AliInfo("Hyperplane fit done\n");
2094
2095                                 // finalize tracklets
2096                                 Int_t labels[12];
2097                                 Int_t outlab[24];
2098                                 Int_t nlab = 0;
2099                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2100                                         if (!cseed[iLayer].IsOK()) continue;
2101
2102                                         if (cseed[iLayer].GetLabels(0) >= 0) {
2103                                                 labels[nlab] = cseed[iLayer].GetLabels(0);
2104                                                 nlab++;
2105                                         }
2106
2107                                         if (cseed[iLayer].GetLabels(1) >= 0) {
2108                                                 labels[nlab] = cseed[iLayer].GetLabels(1);
2109                                                 nlab++;
2110                                         }
2111                                 }
2112                                 Freq(nlab,labels,outlab,kFALSE);
2113                                 Int_t label     = outlab[0];
2114                                 Int_t frequency = outlab[1];
2115                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2116                                         cseed[iLayer].SetFreq(frequency);
2117                                         cseed[iLayer].SetChi2Z(chi2ZF);
2118                                 }
2119
2120                                 if(AliTRDReconstructor::StreamLevel() >= 2){
2121                                         TTreeSRedirector &cstreamer = *fgDebugStreamer;
2122                                         cstreamer << "MakeSeeds2"
2123                                                 << "C="       << curv
2124                                                 << "Chi2TR="  << chi2[0]
2125                                                 << "Chi2TC="  << chi2[1]
2126                                                 << "Chi2RF="  << chi2RF
2127                                                 << "Chi2ZF="  << chi2ZF
2128                                                 << "Nlayers=" << nlayers
2129                                                 << "NUsedS="  << nUsedCl
2130                                                 << "NUsed="   << nusedf
2131                                                 << "Like="    << like
2132                                                 << "S0.="     << &cseed[0]
2133                                                 << "S1.="     << &cseed[1]
2134                                                 << "S2.="     << &cseed[2]
2135                                                 << "S3.="     << &cseed[3]
2136                                                 << "S4.="     << &cseed[4]
2137                                                 << "S5.="     << &cseed[5]
2138                                                 << "Label="   << label
2139                                                 << "Freq="    << frequency
2140                                                 << "\n";
2141                                 }
2142
2143                                 ntracks++;
2144                                 if(ntracks == kMaxTracksStack){
2145                                         AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
2146                                         for(int isl=0; isl<4; isl++) delete layer[isl];
2147                                         return ntracks;
2148                                 }
2149                                 cseed += 6;
2150                         }
2151                 }
2152         }
2153         for(int isl=0; isl<4; isl++) delete layer[isl];
2154
2155         return ntracks;
2156 }
2157
2158 //_____________________________________________________________________________
2159 AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
2160 {
2161   //
2162   // Build a TRD track out of tracklet candidates
2163   //
2164   // Parameters :
2165   //   seeds  : array of tracklets
2166   //   params : track parameters (see MakeSeeds() function body for a detailed description)
2167   //
2168   // Output :
2169   //   The TRD track.
2170   //
2171   // Detailed description
2172   //
2173   // To be discussed with Marian !!
2174   //
2175
2176   Double_t alpha = AliTRDgeometry::GetAlpha();
2177   Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
2178   Double_t c[15];
2179
2180   c[ 0] = 0.2;
2181   c[ 1] = 0.0; c[ 2] = 2.0;
2182   c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
2183   c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0;  c[ 9] = 0.1;
2184   c[10] = 0.0; c[11] = 0.0; c[12] = 0.0;  c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
2185
2186   AliTRDtrackV1 *track = new AliTRDtrackV1(seeds, &params[1], c, params[0], params[6]*alpha+shift);
2187         track->PropagateTo(params[0]-5.0);
2188   track->ResetCovariance(1);
2189   Int_t nc = FollowBackProlongation(*track);
2190         //AliInfo(Form("N clusters for track %d", nc));
2191         if (nc < 30) {
2192     delete track;
2193     track = 0x0;
2194   } else {
2195 //     track->CookdEdx();
2196 //     track->CookdEdxTimBin(-1);
2197 //     CookLabel(track, 0.9);
2198   }
2199
2200   return track;
2201 }
2202
2203 //____________________________________________________________________
2204 void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const
2205 {
2206         // to be implemented, preferably at the level of TRD tracklet. !!!!!!!
2207 }
2208
2209 //____________________________________________________________________
2210 Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDseedV1 *cseed)
2211 {
2212   //
2213   // Sort tracklets according to "quality" and try to "improve" the first 4 worst
2214   //
2215   // Parameters :
2216   //  layers : Array of propagation layers for a stack/supermodule
2217   //  cseed  : Array of 6 seeding tracklets which has to be improved
2218   //
2219   // Output :
2220   //   cssed : Improved seeds
2221   //
2222   // Detailed description
2223   //
2224   // Iterative procedure in which new clusters are searched for each
2225   // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
2226   // can be maximized. If some optimization is found the old seeds are replaced.
2227   //
2228
2229         // make a local working copy
2230         AliTRDtrackingChamber *chamber = 0x0;
2231         AliTRDseedV1 bseed[6];
2232         Int_t nLayers = 0;
2233         for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
2234
2235         Float_t lastquality = 10000.0;
2236         Float_t lastchi2    = 10000.0;
2237         Float_t chi2        =  1000.0;
2238
2239         for (Int_t iter = 0; iter < 4; iter++) {
2240                 Float_t sumquality = 0.0;
2241                 Float_t squality[6];
2242                 Int_t   sortindexes[6];
2243
2244                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
2245                         squality[jLayer]  = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
2246                         sumquality += squality[jLayer];
2247                 }
2248                 if ((sumquality >= lastquality) || (chi2       >     lastchi2)) break;
2249
2250                 nLayers = 0;
2251                 lastquality = sumquality;
2252                 lastchi2    = chi2;
2253                 if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
2254
2255                 TMath::Sort(6, squality, sortindexes, kFALSE);
2256                 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
2257                         Int_t bLayer = sortindexes[jLayer];
2258                         if(!(chamber = stack[bLayer])) continue;
2259                         bseed[bLayer].AttachClustersIter(chamber, squality[bLayer], kTRUE);
2260                         if(bseed[bLayer].IsOK()) nLayers++;
2261                 }
2262
2263                 chi2 = FitTiltedRieman(bseed, kTRUE);
2264         } // Loop: iter
2265
2266         // we are sure that at least 2 tracklets are OK !
2267         return nLayers+2;
2268 }
2269
2270 //_________________________________________________________________________
2271 Double_t AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Double_t *chi2){
2272 //
2273 // Calculates the Track Likelihood value. This parameter serves as main quality criterion for
2274 // the track selection
2275 // The likelihood value containes:
2276 //    - The chi2 values from the both fitters and the chi2 values in z-direction from a linear fit
2277 //    - The Sum of the Parameter  |slope_ref - slope_fit|/Sigma of the tracklets
2278 // For all Parameters an exponential dependency is used
2279 //
2280 // Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate
2281 //             - Array of chi2 values:
2282 //                 * Non-Constrained Tilted Riemann fit
2283 //                 * Vertex-Constrained Tilted Riemann fit
2284 //                 * z-Direction from Linear fit
2285 // Output:     - The calculated track likelihood
2286 //
2287 // debug level 2
2288 //
2289
2290         Double_t sumdaf = 0, nLayers = 0;
2291         for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) {
2292                 if(!tracklets[iLayer].IsOK()) continue;
2293                 sumdaf += TMath::Abs((tracklets[iLayer].GetYfit(1) - tracklets[iLayer].GetYref(1))/ tracklets[iLayer].GetSigmaY2());
2294                 nLayers++;
2295         }
2296         sumdaf /= Float_t (nLayers - 2.0);
2297
2298         Double_t likeChi2Z  = TMath::Exp(-chi2[2] * 0.14);                      // Chi2Z
2299         Double_t likeChi2TC = TMath::Exp(-chi2[1] * 0.677);                     // Constrained Tilted Riemann
2300         Double_t likeChi2TR = TMath::Exp(-chi2[0] * 0.78);                      // Non-constrained Tilted Riemann
2301         Double_t likeAF     = TMath::Exp(-sumdaf * 3.23);
2302         Double_t trackLikelihood     = likeChi2Z * likeChi2TR * likeAF;
2303
2304         if(AliTRDReconstructor::StreamLevel() >= 2){
2305                 TTreeSRedirector &cstreamer = *fgDebugStreamer;
2306                 cstreamer << "CalculateTrackLikelihood0"
2307                         << "LikeChi2Z="         << likeChi2Z
2308                         << "LikeChi2TR="        << likeChi2TR
2309                         << "LikeChi2TC="        << likeChi2TC
2310                         << "LikeAF="                    << likeAF
2311                         << "TrackLikelihood=" << trackLikelihood
2312                         << "\n";
2313         }
2314
2315         return trackLikelihood;
2316 }
2317
2318 //____________________________________________________________________
2319 Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4]
2320                                        , Double_t *chi2)
2321 {
2322   //
2323   // Calculate the probability of this track candidate.
2324   //
2325   // Parameters :
2326   //   cseeds : array of candidate tracklets
2327   //   planes : array of seeding planes (see seeding configuration)
2328   //   chi2   : chi2 values (on the Z and Y direction) from the rieman fit of the track.
2329   //
2330   // Output :
2331   //   likelihood value
2332   //
2333   // Detailed description
2334   //
2335   // The track quality is estimated based on the following 4 criteria:
2336   //  1. precision of the rieman fit on the Y direction (likea)
2337   //  2. chi2 on the Y direction (likechi2y)
2338   //  3. chi2 on the Z direction (likechi2z)
2339   //  4. number of attached clusters compared to a reference value
2340   //     (see AliTRDrecoParam::fkFindable) (likeN)
2341   //
2342   // The distributions for each type of probabilities are given below as of
2343   // (date). They have to be checked to assure consistency of estimation.
2344   //
2345
2346         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2347         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2348         // ratio of the total number of clusters/track which are expected to be found by the tracker.
2349         Float_t fgFindable = AliTRDReconstructor::RecoParam()->GetFindableClusters();
2350
2351
2352         Int_t nclusters = 0;
2353         Double_t sumda = 0.;
2354         for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
2355                 Int_t jlayer = planes[ilayer];
2356                 nclusters += cseed[jlayer].GetN2();
2357                 sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
2358         }
2359         Double_t likea     = TMath::Exp(-sumda*10.6);
2360         Double_t likechi2y  = 0.0000000001;
2361         if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
2362         Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
2363         Int_t enc = Int_t(fgFindable*4.*nTimeBins);     // Expected Number Of Clusters, normally 72
2364         Double_t likeN     = TMath::Exp(-(enc - nclusters) * 0.19);
2365
2366         Double_t like      = likea * likechi2y * likechi2z * likeN;
2367
2368         //AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
2369         if(AliTRDReconstructor::StreamLevel() >= 2){
2370                 TTreeSRedirector &cstreamer = *fgDebugStreamer;
2371                 cstreamer << "CookLikelihood"
2372                         << "sumda="     << sumda
2373                         << "chi0="      << chi2[0]
2374                         << "chi1="      << chi2[1]
2375                         << "likea="     << likea
2376                         << "likechi2y=" << likechi2y
2377                         << "likechi2z=" << likechi2z
2378                         << "nclusters=" << nclusters
2379                         << "likeN="     << likeN
2380                         << "like="      << like
2381                         << "\n";
2382         }
2383
2384         return like;
2385 }
2386
2387
2388 //___________________________________________________________________
2389 void AliTRDtrackerV1::GetMeanCLStack(AliTRDtrackingChamber *chamber, Int_t *planes, Double_t *params)
2390 {
2391   //
2392   // Determines the Mean number of clusters per layer.
2393   // Needed to determine good Seeding Layers
2394   //
2395   // Parameters:
2396   //    - Array of AliTRDchamberTimeBins
2397   //    - Container for the params
2398   //
2399   // Detailed description
2400   //
2401   // Two Iterations:
2402   // In the first Iteration the mean is calculted using all layers.
2403   // After this, all layers outside the 1-sigma-region are rejected.
2404   // Then the mean value and the standard-deviation are calculted a second
2405   // time in order to select all layers in the 1-sigma-region as good-candidates.
2406   //
2407
2408         Float_t mean = 0, stdev = 0;
2409         Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
2410         Int_t position = 0;
2411         memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
2412         memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
2413         Int_t nused = 0;
2414         AliTRDchamberTimeBin *layers = chamber->GetTB(0);
2415         for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
2416                 for(Int_t ils = 0; ils < fTimeBinsPerPlane; ils++){
2417                         position = planes[ipl]*fTimeBinsPerPlane + ils;
2418                         ncl[ipl * fTimeBinsPerPlane + ils] = layers[position].GetNClusters();
2419                         nused = 0;
2420                         for(Int_t icl = 0; icl < ncl[ipl * fTimeBinsPerPlane + ils]; icl++)
2421                                 if((layers[position].GetCluster(icl))->IsUsed()) nused++;
2422                         ncl[ipl * fTimeBinsPerPlane + ils] -= nused;
2423                 }
2424         }
2425         // Declaration of quartils:
2426         //Double_t qvals[3] = {0.0, 0.0, 0.0};
2427         //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
2428         // Iterations
2429         Int_t counter;
2430         Double_t *array;
2431         Int_t *limit;
2432         Int_t nLayers = fTimeBinsPerPlane * kNSeedPlanes;
2433         for(Int_t iter = 0; iter < 2; iter++){
2434                 array = (iter == 0) ? &ncl[0] : &mcl[0];
2435                 limit = (iter == 0) ? &nLayers : &counter;
2436                 counter = 0;
2437                 if(iter == 1){
2438                         for(Int_t i = 0; i < fTimeBinsPerPlane *kNSeedPlanes; i++){
2439                                 if((ncl[i] >  mean + stdev) || (ncl[i] <  mean - stdev)) continue; // Outside 1-sigma region
2440 //                              if((ncl[i] >  qvals[2]) || (ncl[i] <  qvals[0])) continue; // Outside 1-sigma region
2441                                 if(ncl[i] == 0) continue;                                                // 0-Layers also rejected
2442                                 mcl[counter] = ncl[i];
2443                                 counter++;
2444                         }
2445                 }
2446                 if(*limit == 0) break;
2447                 printf("Limit = %d\n", *limit);
2448                 //using quartils instead of mean and RMS
2449 //              TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
2450                 mean = TMath::Median(*limit, array, 0x0);
2451                 stdev  = TMath::RMS(*limit, array);
2452         }
2453 //      printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
2454 //      memcpy(params,qvals,sizeof(Double_t)*3);
2455         params[1] = (Double_t)TMath::Nint(mean);
2456         params[0] = (Double_t)TMath::Nint(mean - stdev);
2457         params[2] = (Double_t)TMath::Nint(mean + stdev);
2458
2459 }
2460
2461 //___________________________________________________________________
2462 Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDtrackingChamber *chamber, Double_t *params)
2463 {
2464   //
2465   // Algorithm to find optimal seeding layer
2466   // Layers inside one sigma region (given by Quantiles) are sorted
2467   // according to their difference.
2468   // All layers outside are sorted according t
2469   //
2470   // Parameters:
2471   //     - Array of AliTRDchamberTimeBins (in the current plane !!!)
2472   //     - Container for the Indices of the seeding Layer candidates
2473   //
2474   // Output:
2475   //     - Number of Layers inside the 1-sigma-region
2476   //
2477   // The optimal seeding layer should contain the mean number of
2478   // custers in the layers in one chamber.
2479   //
2480
2481         //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
2482         const Int_t kMaxClustersLayer = AliTRDchamberTimeBin::kMaxClustersLayer;
2483         Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
2484         memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
2485         memset(indices, 0, sizeof(Int_t)*kNTimeBins);
2486         memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
2487
2488         AliTRDchamberTimeBin *layers = chamber->GetTB(0);
2489         Int_t nused = 0;
2490         for(Int_t ils = 0; ils < fTimeBinsPerPlane; ils++){
2491                 ncl[ils] = layers[ils].GetNClusters();
2492                 nused = 0;
2493                 for(Int_t icl = 0; icl < ncl[ils]; icl++)
2494                         if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
2495                 ncl[ils] -= nused;
2496         }
2497
2498         Float_t mean = params[1];
2499         for(Int_t ils = 0; ils < fTimeBinsPerPlane; ils++){
2500                 memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(fTimeBinsPerPlane - ils));
2501                 indices[bins[ncl[ils]+1]] = ils;
2502                 for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
2503                         bins[i]++;
2504         }
2505
2506         //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
2507         Int_t sbin = -1;
2508         Int_t nElements;
2509         Int_t position = 0;
2510         TRandom *r = new TRandom();
2511         Int_t iter = 0;
2512         while(1){
2513                 while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
2514                         // Randomly selecting one bin
2515                         sbin = (Int_t)r->Poisson(mean);
2516                 }
2517                 printf("Bin = %d\n",sbin);
2518                 //Randomly selecting one Layer in the bin
2519                 nElements = bins[sbin + 1] - bins[sbin];
2520                 printf("nElements = %d\n", nElements);
2521                 if(iter == 5){
2522                         position = (Int_t)(gRandom->Rndm()*(fTimeBinsPerPlane-1));
2523                         break;
2524                 }
2525                 else if(nElements==0){
2526                         iter++;
2527                         continue;
2528                 }
2529                 position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
2530                 break;
2531         }
2532         delete r;
2533         return indices[position];
2534 }
2535
2536 //____________________________________________________________________
2537 AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDtrackingChamber *chamber, AliTRDseedV1* reference) const
2538 {
2539   //
2540   // Finds a seeding Cluster for the extrapolation chamber.
2541   //
2542   // The seeding cluster should be as close as possible to the assumed
2543   // track which is represented by a Rieman fit.
2544   // Therefore the selecting criterion is the minimum distance between
2545   // the best fitting cluster and the Reference which is derived from
2546   // the AliTRDseed. Because all layers are assumed to be equally good
2547   // a linear search is performed.
2548   //
2549   // Imput parameters: - layers: array of AliTRDchamberTimeBins (in one chamber!!!)
2550   //                   - sfit: the reference
2551   //
2552   // Output:           - the best seeding cluster
2553   //
2554
2555
2556         // distances as squared distances
2557         Int_t index = 0;
2558         Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearestDistance =100000.0;
2559         ypos = reference->GetYref(0);
2560         zpos = reference->GetZref(0);
2561         AliTRDcluster *currentBest = 0x0, *temp = 0x0;
2562         AliTRDchamberTimeBin *layers = chamber->GetTB(0);
2563         for(Int_t ils = 0; ils < fTimeBinsPerPlane; ils++){
2564                 // Reference positions
2565 //              ypos = reference->GetYat(layers[ils].GetX());
2566 //              zpos = reference->GetZat(layers[ils].GetX());
2567                 index = layers[ils].SearchNearestCluster(ypos, zpos, AliTRDReconstructor::RecoParam()->GetRoad2y(), AliTRDReconstructor::RecoParam()->GetRoad2z());
2568                 if(index == -1) continue;
2569                 temp = layers[ils].GetCluster(index);
2570                 if(!temp) continue;
2571                 distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
2572                 if(distance < nearestDistance){
2573                         nearestDistance = distance;
2574                         currentBest = temp;
2575                 }
2576         }
2577         return currentBest;
2578 }
2579
2580
2581 //____________________________________________________________________
2582 void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4])
2583 {
2584   //
2585   // Map seeding configurations to detector planes.
2586   //
2587   // Parameters :
2588   //   iconfig : configuration index
2589   //   planes  : member planes of this configuration. On input empty.
2590   //
2591   // Output :
2592   //   planes : contains the planes which are defining the configuration
2593   //
2594   // Detailed description
2595   //
2596   // Here is the list of seeding planes configurations together with
2597   // their topological classification:
2598   //
2599   //  0 - 5432 TQ 0
2600   //  1 - 4321 TQ 0
2601   //  2 - 3210 TQ 0
2602   //  3 - 5321 TQ 1
2603   //  4 - 4210 TQ 1
2604   //  5 - 5431 TQ 1
2605   //  6 - 4320 TQ 1
2606   //  7 - 5430 TQ 2
2607   //  8 - 5210 TQ 2
2608   //  9 - 5421 TQ 3
2609   // 10 - 4310 TQ 3
2610   // 11 - 5410 TQ 4
2611   // 12 - 5420 TQ 5
2612   // 13 - 5320 TQ 5
2613   // 14 - 5310 TQ 5
2614   //
2615   // The topologic quality is modeled as follows:
2616   // 1. The general model is define by the equation:
2617   //  p(conf) = exp(-conf/2)
2618   // 2. According to the topologic classification, configurations from the same
2619   //    class are assigned the agerage value over the model values.
2620   // 3. Quality values are normalized.
2621   //
2622   // The topologic quality distribution as function of configuration is given below:
2623   //Begin_Html
2624   // <img src="gif/topologicQA.gif">
2625   //End_Html
2626   //
2627
2628         switch(iconfig){
2629         case 0: // 5432 TQ 0
2630                 planes[0] = 2;
2631                 planes[1] = 3;
2632                 planes[2] = 4;
2633                 planes[3] = 5;
2634                 break;
2635         case 1: // 4321 TQ 0
2636                 planes[0] = 1;
2637                 planes[1] = 2;
2638                 planes[2] = 3;
2639                 planes[3] = 4;
2640                 break;
2641         case 2: // 3210 TQ 0
2642                 planes[0] = 0;
2643                 planes[1] = 1;
2644                 planes[2] = 2;
2645                 planes[3] = 3;
2646                 break;
2647         case 3: // 5321 TQ 1
2648                 planes[0] = 1;
2649                 planes[1] = 2;
2650                 planes[2] = 3;
2651                 planes[3] = 5;
2652                 break;
2653         case 4: // 4210 TQ 1
2654                 planes[0] = 0;
2655                 planes[1] = 1;
2656                 planes[2] = 2;
2657                 planes[3] = 4;
2658                 break;
2659         case 5: // 5431 TQ 1
2660                 planes[0] = 1;
2661                 planes[1] = 3;
2662                 planes[2] = 4;
2663                 planes[3] = 5;
2664                 break;
2665         case 6: // 4320 TQ 1
2666                 planes[0] = 0;
2667                 planes[1] = 2;
2668                 planes[2] = 3;
2669                 planes[3] = 4;
2670                 break;
2671         case 7: // 5430 TQ 2
2672                 planes[0] = 0;
2673                 planes[1] = 3;
2674                 planes[2] = 4;
2675                 planes[3] = 5;
2676                 break;
2677         case 8: // 5210 TQ 2
2678                 planes[0] = 0;
2679                 planes[1] = 1;
2680                 planes[2] = 2;
2681                 planes[3] = 5;
2682                 break;
2683         case 9: // 5421 TQ 3
2684                 planes[0] = 1;
2685                 planes[1] = 2;
2686                 planes[2] = 4;
2687                 planes[3] = 5;
2688                 break;
2689         case 10: // 4310 TQ 3
2690                 planes[0] = 0;
2691                 planes[1] = 1;
2692                 planes[2] = 3;
2693                 planes[3] = 4;
2694                 break;
2695         case 11: // 5410 TQ 4
2696                 planes[0] = 0;
2697                 planes[1] = 1;
2698                 planes[2] = 4;
2699                 planes[3] = 5;
2700                 break;
2701         case 12: // 5420 TQ 5
2702                 planes[0] = 0;
2703                 planes[1] = 2;
2704                 planes[2] = 4;
2705                 planes[3] = 5;
2706                 break;
2707         case 13: // 5320 TQ 5
2708                 planes[0] = 0;
2709                 planes[1] = 2;
2710                 planes[2] = 3;
2711                 planes[3] = 5;
2712                 break;
2713         case 14: // 5310 TQ 5
2714                 planes[0] = 0;
2715                 planes[1] = 1;
2716                 planes[2] = 3;
2717                 planes[3] = 5;
2718                 break;
2719         }
2720 }
2721
2722 //____________________________________________________________________
2723 void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
2724 {
2725   //
2726   // Returns the extrapolation planes for a seeding configuration.
2727   //
2728   // Parameters :
2729   //   iconfig : configuration index
2730   //   planes  : planes which are not in this configuration. On input empty.
2731   //
2732   // Output :
2733   //   planes : contains the planes which are not in the configuration
2734   //
2735   // Detailed description
2736   //
2737
2738         switch(iconfig){
2739         case 0: // 5432 TQ 0
2740                 planes[0] = 1;
2741                 planes[1] = 0;
2742                 break;
2743         case 1: // 4321 TQ 0
2744                 planes[0] = 5;
2745                 planes[1] = 0;
2746                 break;
2747         case 2: // 3210 TQ 0
2748                 planes[0] = 4;
2749                 planes[1] = 5;
2750                 break;
2751         case 3: // 5321 TQ 1
2752                 planes[0] = 4;
2753                 planes[1] = 0;
2754                 break;
2755         case 4: // 4210 TQ 1
2756                 planes[0] = 5;
2757                 planes[1] = 3;
2758                 break;
2759         case 5: // 5431 TQ 1
2760                 planes[0] = 2;
2761                 planes[1] = 0;
2762                 break;
2763         case 6: // 4320 TQ 1
2764                 planes[0] = 5;
2765                 planes[1] = 1;
2766                 break;
2767         case 7: // 5430 TQ 2
2768                 planes[0] = 2;
2769                 planes[1] = 1;
2770                 break;
2771         case 8: // 5210 TQ 2
2772                 planes[0] = 4;
2773                 planes[1] = 3;
2774                 break;
2775         case 9: // 5421 TQ 3
2776                 planes[0] = 3;
2777                 planes[1] = 0;
2778                 break;
2779         case 10: // 4310 TQ 3
2780                 planes[0] = 5;
2781                 planes[1] = 2;
2782                 break;
2783         case 11: // 5410 TQ 4
2784                 planes[0] = 3;
2785                 planes[1] = 2;
2786                 break;
2787         case 12: // 5420 TQ 5
2788                 planes[0] = 3;
2789                 planes[1] = 1;
2790                 break;
2791         case 13: // 5320 TQ 5
2792                 planes[0] = 4;
2793                 planes[1] = 1;
2794                 break;
2795         case 14: // 5310 TQ 5
2796                 planes[0] = 4;
2797                 planes[1] = 2;
2798                 break;
2799         }
2800 }
2801
2802 //____________________________________________________________________
2803 AliCluster* AliTRDtrackerV1::GetCluster(Int_t idx) const
2804 {
2805         Int_t ncls = fClusters->GetEntriesFast();
2806         return idx >= 0 || idx < ncls ? (AliCluster*)fClusters->UncheckedAt(idx) : 0x0;
2807 }
2808
2809
2810 //_____________________________________________________________________________
2811 Int_t AliTRDtrackerV1::Freq(Int_t n, const Int_t *inlist
2812                         , Int_t *outlist, Bool_t down)
2813 {
2814   //
2815   // Sort eleements according occurancy
2816   // The size of output array has is 2*n
2817   //
2818
2819   if (n <= 0) {
2820     return 0;
2821   }
2822
2823   Int_t *sindexS = new Int_t[n];   // Temporary array for sorting
2824   Int_t *sindexF = new Int_t[2*n];
2825   for (Int_t i = 0; i < n; i++) {
2826     sindexF[i] = 0;
2827   }
2828
2829   TMath::Sort(n,inlist,sindexS,down);
2830
2831   Int_t last     = inlist[sindexS[0]];
2832   Int_t val      = last;
2833   sindexF[0]     = 1;
2834   sindexF[0+n]   = last;
2835   Int_t countPos = 0;
2836
2837   // Find frequency
2838   for (Int_t i = 1; i < n; i++) {
2839     val = inlist[sindexS[i]];
2840     if (last == val) {
2841       sindexF[countPos]++;
2842     }
2843     else {
2844       countPos++;
2845       sindexF[countPos+n] = val;
2846       sindexF[countPos]++;
2847       last                = val;
2848     }
2849   }
2850   if (last == val) {
2851     countPos++;
2852   }
2853
2854   // Sort according frequency
2855   TMath::Sort(countPos,sindexF,sindexS,kTRUE);
2856
2857   for (Int_t i = 0; i < countPos; i++) {
2858     outlist[2*i  ] = sindexF[sindexS[i]+n];
2859     outlist[2*i+1] = sindexF[sindexS[i]];
2860   }
2861
2862   delete [] sindexS;
2863   delete [] sindexF;
2864
2865   return countPos;
2866
2867 }