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 <TObjArray.h>
  40 #include <TROOT.h>
  41 #include <TTree.h>
  42 #include <TClonesArray.h>
  43 #include <TRandom.h>
  44 #include <TTreeStream.h>
  45
  46 #include "AliLog.h"
  47 #include "AliESDEvent.h"
  48 #include "AliAlignObj.h"
  49 #include "AliRieman.h"
  50 #include "AliTrackPointArray.h"
  51
  52 #include "AliTRDtracker.h"
  53 #include "AliTRDtrackerV1.h"
  54 #include "AliTRDgeometry.h"
  55 #include "AliTRDpadPlane.h"
  56 #include "AliTRDgeometry.h"
  57 #include "AliTRDcluster.h"
  58 #include "AliTRDtrack.h"
  59 #include "AliTRDseed.h"
  60 #include "AliTRDcalibDB.h"
  61 #include "AliTRDCommonParam.h"
  62 #include "AliTRDReconstructor.h"
  63 #include "AliTRDCalibraFillHisto.h"
  64 #include "AliTRDtrackerFitter.h"
  65 #include "AliTRDstackLayer.h"
  66 #include "AliTRDrecoParam.h"
  67 #include "AliTRDseedV1.h"
  68 #include "AliTRDtrackV1.h"
  69 #include "Cal/AliTRDCalDet.h"
  70
  71 #define DEBUG
  72
  73 ClassImp(AliTRDtrackerV1)
  74 Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
  75                 0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
  76                 0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
  77                 0.0474, 0.0408, 0.0335, 0.0335, 0.0335
  78 };
  79
  80 //____________________________________________________________________
  81 AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p)
  82   :AliTRDtracker()
  83   ,fSieveSeeding(0)
  84   ,fTracklets(0x0)
  85   ,fRecoParam(p)
  86   ,fFitter(0x0)
  87 {
  88   //
  89   // Default constructor. Nothing is initialized.
  90   //
  91
  92 }
  93
  94 //____________________________________________________________________
  95 AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p)
  96   :AliTRDtracker(in)
  97   ,fSieveSeeding(0)
  98   ,fTracklets(0x0)
  99   ,fRecoParam(p)
 100   ,fFitter(0x0)
 101 {
 102   //
 103   // Standard constructor.
 104   // Setting of the geometry file, debug output (if enabled)
 105   // and initilize fitter helper.
 106   //
 107
 108         fFitter = new AliTRDtrackerFitter();
 109
 110 #ifdef DEBUG
 111         fFitter->SetDebugStream(fDebugStreamer);
 112 #endif
 113
 114 }
 115
 116 //____________________________________________________________________
 117 AliTRDtrackerV1::~AliTRDtrackerV1()
 118 {
 119   //
 120   // Destructor
 121   //
 122
 123         if(fFitter) delete fFitter;
 124         if(fRecoParam) delete fRecoParam;
 125         if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
 126 }
 127
 128 //____________________________________________________________________
 129 Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
 130 {
 131   //
 132   // Steering stand alone tracking for full TRD detector
 133   //
 134   // Parameters :
 135   //   esd     : The ESD event. On output it contains
 136   //             the ESD tracks found in TRD.
 137   //
 138   // Output :
 139   //   Number of tracks found in the TRD detector.
 140   //
 141   // Detailed description
 142   // 1. Launch individual SM trackers.
 143   //    See AliTRDtrackerV1::Clusters2TracksSM() for details.
 144   //
 145
 146         if(!fRecoParam){
 147                 AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
 148                 return 0;
 149         }
 150
 151         //AliInfo("Start Track Finder ...");
 152         Int_t ntracks = 0;
 153         for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
 154                         //AliInfo(Form("Processing supermodule %i ...", ism));
 155                         ntracks += Clusters2TracksSM(fTrSec[ism], esd);
 156         }
 157         AliInfo(Form("Found %d TRD tracks.", ntracks));
 158         return ntracks;
 159 }
 160
 161
 162 //_____________________________________________________________________________
 163 Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t /*index*/, AliTrackPoint &/*p*/) const
 164 {
 165         //AliInfo(Form("Asking for tracklet %d", index));
 166
 167         if(index<0) return kFALSE;
 168         //AliTRDseedV1 *tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index);
 169         // etc
 170         return kTRUE;
 171 }
 172
 173
 174 //_____________________________________________________________________________
 175 Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
 176 {
 177   //
 178   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
 179   // backpropagated by the TPC tracker. Each seed is first propagated
 180   // to the TRD, and then its prolongation is searched in the TRD.
 181   // If sufficiently long continuation of the track is found in the TRD
 182   // the track is updated, otherwise it's stored as originaly defined
 183   // by the TPC tracker.
 184   //
 185
 186         Int_t   found    = 0;     // number of tracks found
 187         Float_t foundMin = 20.0;
 188
 189         AliTRDseed::SetNTimeBins(fTimeBinsPerPlane);
 190         Int_t    nSeed   = event->GetNumberOfTracks();
 191         if(!nSeed){
 192                 // run stand alone tracking
 193                 if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
 194                 return 0;
 195         }
 196
 197         Float_t *quality = new Float_t[nSeed];
 198         Int_t   *index   = new Int_t[nSeed];
 199         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
 200                 AliESDtrack *seed = event->GetTrack(iSeed);
 201                 Double_t covariance[15];
 202                 seed->GetExternalCovariance(covariance);
 203                 quality[iSeed] = covariance[0] + covariance[2];
 204         }
 205         // Sort tracks according to covariance of local Y and Z
 206         TMath::Sort(nSeed,quality,index,kFALSE);
 207
 208         // Backpropagate all seeds
 209         for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
 210
 211                 // Get the seeds in sorted sequence
 212                 AliESDtrack *seed = event->GetTrack(index[iSeed]);
 213
 214                 // Check the seed status
 215                 ULong_t status = seed->GetStatus();
 216                 if ((status & AliESDtrack::kTPCout) == 0) continue;
 217                 if ((status & AliESDtrack::kTRDout) != 0) continue;
 218
 219                 // Do the back prolongation
 220                 Int_t   lbl         = seed->GetLabel();
 221                 AliTRDtrackV1 *track  = new AliTRDtrackV1(*seed);
 222                 //track->Print();
 223                 track->SetSeedLabel(lbl);
 224                 seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); // Make backup
 225                 fNseeds++;
 226                 Float_t p4          = track->GetC();
 227                 Int_t   expectedClr = FollowBackProlongation(*track);
 228                 //AliInfo(Form("\nTRACK %d Clusters %d [%d] in chi2 %f", index[iSeed], expectedClr, track->GetNumberOfClusters(), track->GetChi2()));
 229                 //track->Print();
 230
 231                 //Double_t cov[15];
 232                 //seed->GetExternalCovariance(cov);
 233                 //AliInfo(Form("track %d cov[%f %f] 0", index[iSeed], cov[0], cov[2]));
 234
 235                 if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
 236                                 (track->Pt() > 0.8)) {
 237                         //
 238                         // Make backup for back propagation
 239                         //
 240                         Int_t foundClr = track->GetNumberOfClusters();
 241                         if (foundClr >= foundMin) {
 242                                 //AliInfo(Form("Making backup track ncls [%d]...", foundClr));
 243                                 track->CookdEdx();
 244                                 track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
 245                                 CookLabel(track,1 - fgkLabelFraction);
 246                                 if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
 247
 248
 249                 //seed->GetExternalCovariance(cov);
 250                 //AliInfo(Form("track %d cov[%f %f] 0 test", index[iSeed], cov[0], cov[2]));
 251
 252                                 // Sign only gold tracks
 253                                 if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
 254                                         if ((seed->GetKinkIndex(0)      ==   0) &&
 255                                                         (track->Pt()                <  1.5)) UseClusters(track);
 256                                 }
 257                                 Bool_t isGold = kFALSE;
 258
 259                                 // Full gold track
 260                                 if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
 261                                         if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 262
 263                                         isGold = kTRUE;
 264                                 }
 265                 //seed->GetExternalCovariance(cov);
 266                 //AliInfo(Form("track %d cov[%f %f] 00", index[iSeed], cov[0], cov[2]));
 267
 268                                 // Almost gold track
 269                                 if ((!isGold)  && (track->GetNCross() == 0) &&
 270                                                 (track->GetChi2() / track->GetNumberOfClusters()  < 7)) {
 271                                         //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 272                                         if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 273
 274                                         isGold = kTRUE;
 275                                 }
 276                 //seed->GetExternalCovariance(cov);
 277                 //AliInfo(Form("track %d cov[%f %f] 01", index[iSeed], cov[0], cov[2]));
 278
 279                                 if ((!isGold) && (track->GetBackupTrack())) {
 280                                         if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
 281                                                 seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
 282                                                 isGold = kTRUE;
 283                                         }
 284                                 }
 285                 //seed->GetExternalCovariance(cov);
 286                 //AliInfo(Form("track %d cov[%f %f] 02", index[iSeed], cov[0], cov[2]));
 287
 288                                 //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected())  > 0.4)) {
 289                                         //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 290                                 //}
 291                         }
 292                 }
 293                 /**/
 294
 295                 /**/
 296                 // Debug part of tracking
 297 /*              TTreeSRedirector &cstream = *fDebugStreamer;
 298                 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.
 299                 if (AliTRDReconstructor::StreamLevel() > 0) {
 300                         if (track->GetBackupTrack()) {
 301                                 cstream << "Tracks"
 302                                 << "EventNrInFile="  << eventNrInFile
 303                                 << "ESD.="     << seed
 304                                 << "trd.="     << track
 305                                 << "trdback.=" << track->GetBackupTrack()
 306                                 << "\n";
 307                         }
 308                         else {
 309                                 cstream << "Tracks"
 310                                 << "EventNrInFile="  << eventNrInFile
 311                                 << "ESD.="     << seed
 312                                 << "trd.="     << track
 313                                 << "trdback.=" << track
 314                                 << "\n";
 315                         }
 316                 }*/
 317                 /**/
 318
 319                 //seed->GetExternalCovariance(cov);
 320                 //AliInfo(Form("track %d cov[%f %f] 1", index[iSeed], cov[0], cov[2]));
 321
 322                 // Propagation to the TOF (I.Belikov)
 323                 if (track->GetStop() == kFALSE) {
 324                         //AliInfo("Track not stopped in TRD ...");
 325                         Double_t xtof  = 371.0;
 326                         Double_t xTOF0 = 370.0;
 327
 328                         Double_t c2    = track->GetSnp() + track->GetC() * (xtof - track->GetX());
 329                         if (TMath::Abs(c2) >= 0.99) {
 330                                 delete track;
 331                                 continue;
 332                         }
 333
 334                         PropagateToX(*track,xTOF0,fgkMaxStep);
 335
 336                         // Energy losses taken to the account - check one more time
 337                         c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
 338                         if (TMath::Abs(c2) >= 0.99) {
 339                                 delete track;
 340                                 continue;
 341                         }
 342
 343                         //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
 344                         //      fHBackfit->Fill(7);
 345                         //delete track;
 346                         //      continue;
 347                         //}
 348
 349                         Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
 350                         Double_t y;
 351                         track->GetYAt(xtof,GetBz(),y);
 352                         if (y >  ymax) {
 353                                 if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
 354                                         delete track;
 355                                         continue;
 356                                 }
 357                         }else if (y < -ymax) {
 358                                 if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
 359                                         delete track;
 360                                         continue;
 361                                 }
 362                         }
 363
 364                         if (track->PropagateTo(xtof)) {
 365                                 //AliInfo("set kTRDout");
 366                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
 367
 368                                 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 369                                         for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
 370                                                 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
 371                                         }
 372                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 373                                 }
 374                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
 375                                 if (track->GetNumberOfClusters() > foundMin) found++;
 376                         }
 377                 } else {
 378                         //AliInfo("Track stopped in TRD ...");
 379
 380                         if ((track->GetNumberOfClusters() >              15) &&
 381                                         (track->GetNumberOfClusters() > 0.5*expectedClr)) {
 382                                 seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
 383
 384                                 //seed->SetStatus(AliESDtrack::kTRDStop);
 385                                 for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
 386                                         for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
 387                                                 seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
 388                                         }
 389                                         seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 390                                 }
 391                                 //seed->SetTRDtrack(new AliTRDtrack(*track));
 392                                 found++;
 393                         }
 394                 }
 395
 396                 //if (((t->GetStatus()&AliESDtrack::kTRDout)!=0 )
 397
 398                 seed->SetTRDQuality(track->StatusForTOF());
 399                 seed->SetTRDBudget(track->GetBudget(0));
 400
 401 //              if ((seed->GetStatus()&AliESDtrack::kTRDin)!=0 ) printf("TRDin ");
 402 //              if ((seed->GetStatus()&AliESDtrack::kTRDbackup)!=0 ) printf("TRDbackup ");
 403 //              if ((seed->GetStatus()&AliESDtrack::kTRDStop)!=0 ) printf("TRDstop ");
 404 //              if ((seed->GetStatus()&AliESDtrack::kTRDout)!=0 ) printf("TRDout ");
 405 //              printf("\n");
 406                 delete track;
 407
 408                 //seed->GetExternalCovariance(cov);
 409                 //AliInfo(Form("track %d cov[%f %f] 2", index[iSeed], cov[0], cov[2]));
 410         }
 411
 412
 413         AliInfo(Form("Number of seeds: %d",fNseeds));
 414         AliInfo(Form("Number of back propagated TRD tracks: %d",found));
 415
 416         //fSeeds->Clear();
 417         fNseeds = 0;
 418
 419         delete [] index;
 420         delete [] quality;
 421
 422   return 0;
 423 }
 424
 425
 426 //____________________________________________________________________
 427 Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event)
 428 {
 429   //
 430   // Refits tracks within the TRD. The ESD event is expected to contain seeds
 431   // at the outer part of the TRD.
 432   // The tracks are propagated to the innermost time bin
 433   // of the TRD and the ESD event is updated
 434   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
 435   //
 436
 437   Int_t   nseed    = 0; // contor for loaded seeds
 438   Int_t   found    = 0; // contor for updated TRD tracks
 439
 440   // Calibration monitor
 441   AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
 442   if (!calibra) AliInfo("Could not get Calibra instance\n");
 443
 444
 445   AliTRDtrackV1 track;
 446   for (Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++) {
 447     AliESDtrack *seed = event->GetTrack(itrack);
 448                 new(&track) AliTRDtrackV1(*seed);
 449
 450     if (track.GetX() < 270.0) {
 451       seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup);
 452       //AliInfo(Form("Remove for X = %7.3f [270.]\n", track.GetX()));
 453                         continue;
 454     }
 455
 456     ULong_t status = seed->GetStatus();
 457     if((status & AliESDtrack::kTRDout) == 0) continue;
 458     if((status & AliESDtrack::kTRDin)  != 0) continue;
 459     nseed++;
 460
 461     track.ResetCovariance(50.0);
 462
 463                 // do the propagation and processing
 464     Bool_t kUPDATE = kFALSE;
 465                 Double_t xTPC = 250.0;
 466     if(FollowProlongation(track)){
 467                         // computes PID for track
 468                         track.CookPID();
 469                         // update calibration references using this track
 470                         if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(&track);
 471
 472                         // Prolongate to TPC
 473                         if (PropagateToX(track, xTPC, fgkMaxStep)) { //  -with update
 474                                 seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit);
 475                                 track.UpdateESDtrack(seed);
 476                                 // Add TRD track to ESDfriendTrack
 477                                 if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
 478                                         AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
 479                                         calibTrack->SetOwner();
 480                                         seed->AddCalibObject(calibTrack);
 481                                 }
 482                                 found++;
 483                                 kUPDATE = kTRUE;
 484                         }
 485                 }
 486
 487                 // Prolongate to TPC without update
 488                 if(!kUPDATE) {
 489       AliTRDtrackV1 tt(*seed);
 490       if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDrefit);
 491     }
 492   }
 493   AliInfo(Form("Number of loaded seeds: %d",nseed));
 494   AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
 495
 496         return 0;
 497 }
 498
 499
 500 //____________________________________________________________________
 501 Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
 502 {
 503 // Extrapolates the TRD track in the TPC direction.
 504 //
 505 // Parameters
 506 //   t : the TRD track which has to be extrapolated
 507 //
 508 // Output
 509 //   number of clusters attached to the track
 510 //
 511 // Detailed description
 512 //
 513 // Starting from current radial position of track <t> this function
 514 // extrapolates the track through the 6 TRD layers. The following steps
 515 // are being performed for each plane:
 516 // 1. prepare track:
 517 //   a. get plane limits in the local x direction
 518 //   b. check crossing sectors
 519 //   c. check track inclination
 520 // 2. search tracklet in the tracker list (see GetTracklet() for details)
 521 // 3. evaluate material budget using the geo manager
 522 // 4. propagate and update track using the tracklet information.
 523 //
 524 // Debug level 2
 525 //
 526
 527         //AliInfo("");
 528         Int_t    nClustersExpected = 0;
 529         Float_t clength = AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick();
 530         Int_t lastplane = 5; //GetLastPlane(&t);
 531         for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
 532                 //AliInfo(Form("plane %d", iplane));
 533     Int_t row1 = GetGlobalTimeBin(0, iplane, 0); // to be modified to the true time bin in the geometrical acceptance
 534                 //AliInfo(Form("row1 %d", row1));
 535
 536     // Propagate track close to the plane if neccessary
 537                 AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row1);
 538                 Double_t currentx  = layer->GetX();
 539     if (currentx < (-fgkMaxStep + t.GetX()))
 540       if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
 541
 542     if (!AdjustSector(&t)) break;
 543
 544                 Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
 545                 //AliInfo(Form("row0 %d", row0));
 546
 547     // Start global position
 548     Double_t xyz0[3];
 549     t.GetXYZ(xyz0);
 550
 551                 // End global position
 552     Double_t x = fTrSec[0]->GetLayer(row0)->GetX(), y, z;
 553     if (!t.GetProlongation(x,y,z)) break;
 554     Double_t xyz1[3];
 555     xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
 556     xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
 557     xyz1[2] =  z;
 558
 559     // Get material budget
 560     Double_t param[7];
 561     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
 562     Double_t xrho= param[0]*param[4];
 563     Double_t xx0 = param[1]; // Get mean propagation parameters
 564
 565     // Propagate and update
 566     //Int_t sector = t.GetSector();
 567     Int_t   index   = 0;
 568                 //AliInfo(Form("sector %d", sector));
 569     AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
 570           //AliInfo(Form("tracklet %p @ %d", tracklet, index));
 571                 if(!tracklet) continue;
 572                 //AliInfo(Form("reco %p", tracklet->GetRecoParam()));
 573                 t.SetTracklet(tracklet, iplane, index);
 574
 575                 t.PropagateTo(tracklet->GetX0() - clength, xx0, xrho);
 576           if (!AdjustSector(&t)) break;
 577
 578     Double_t maxChi2 = t.GetPredictedChi2(tracklet);
 579           if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){
 580                 nClustersExpected += tracklet->GetN();
 581         }
 582   }
 583
 584 #ifdef DEBUG
 585         if(AliTRDReconstructor::StreamLevel() > 1){
 586                 Int_t index;
 587                 for(int iplane=0; iplane<6; iplane++){
 588                         AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
 589                         if(!tracklet) continue;
 590                         t.SetTracklet(tracklet, iplane, index);
 591                 }
 592
 593                 TTreeSRedirector &cstreamer = *fDebugStreamer;
 594                 cstreamer << "FollowProlongation"
 595                         << "ncl="      << nClustersExpected
 596                         << "track.="   << &t
 597                         << "\n";
 598         }
 599 #endif
 600
 601   return nClustersExpected;
 602
 603 }
 604
 605 //_____________________________________________________________________________
 606 Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
 607 {
 608 // Extrapolates the TRD track in the TOF direction.
 609 //
 610 // Parameters
 611 //   t : the TRD track which has to be extrapolated
 612 //
 613 // Output
 614 //   number of clusters attached to the track
 615 //
 616 // Detailed description
 617 //
 618 // Starting from current radial position of track <t> this function
 619 // extrapolates the track through the 6 TRD layers. The following steps
 620 // are being performed for each plane:
 621 // 1. prepare track:
 622 //   a. get plane limits in the local x direction
 623 //   b. check crossing sectors
 624 //   c. check track inclination
 625 // 2. build tracklet (see AliTRDseed::AttachClusters() for details)
 626 // 3. evaluate material budget using the geo manager
 627 // 4. propagate and update track using the tracklet information.
 628 //
 629 // Debug level 2
 630 //
 631
 632         Int_t nClustersExpected = 0;
 633   // Loop through the TRD planes
 634   for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
 635                 //AliInfo(Form("Processing plane %d ...", iplane));
 636                 // Get the global time bin for the first local time bin of the given plane
 637     Int_t row0 = GetGlobalTimeBin(0, iplane, fTimeBinsPerPlane-1);
 638
 639                 // Retrive first propagation layer in the chamber
 640                 AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row0);
 641
 642                 // Get the X coordinates of the propagation layer for the first time bin
 643     Double_t currentx = layer->GetX(); // what if X is not defined ???
 644     if (currentx < t.GetX()) continue;
 645
 646                 // Get the global time bin for the last local time bin of the given plane
 647                 Int_t row1 = GetGlobalTimeBin(0, iplane, 0);
 648
 649                 // Propagate closer to the current chamber if neccessary
 650     if (currentx > (fgkMaxStep + t.GetX()) && !PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
 651
 652                 // Rotate track to adjacent sector if neccessary
 653     if (!AdjustSector(&t)) break;
 654                 Int_t sector =  Int_t(TMath::Abs(t.GetAlpha()/AliTRDgeometry::GetAlpha()));
 655                 if(t.GetAlpha() < 0) sector = AliTRDgeometry::Nsect() - sector-1;
 656
 657                 //AliInfo(Form("sector %d [%f]", sector, t.GetAlpha()));
 658
 659                 // Check whether azimuthal angle is getting too large
 660     if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
 661
 662     //Calculate global entry and exit positions of the track in chamber (only track prolongation)
 663     Double_t xyz0[3]; // entry point
 664                 t.GetXYZ(xyz0);
 665                 //printf("Entry global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz0[0], xyz0[1], xyz0[2]);
 666
 667                 // Get local Y and Z at the X-position of the end of the chamber
 668                 Double_t x0 = fTrSec[sector]->GetLayer(row1)->GetX(), y, z;
 669                 if (!t.GetProlongation(x0, y, z)) break;
 670                 //printf("Exit  local x[%7.3f] y[%7.3f] z[%7.3f]\n", x0, y, z);
 671
 672                 Double_t xyz1[3]; // exit point
 673                 xyz1[0] =  x0 * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
 674     xyz1[1] = +x0 * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
 675     xyz1[2] =  z;
 676
 677                 //printf("Exit  global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz1[0], xyz1[1], xyz1[2]);
 678                 // Find tracklet along the path inside the chamber
 679     AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
 680                 // if the track is not already build (e.g. stand alone tracker) we build it now.
 681                 if(!tracklet.GetN()){ // a better check has to be implemented TODO!!!!!!!
 682
 683                         //AliInfo(Form("Building tracklet for plane %d ...", iplane));
 684                         // check if we are inside detection volume
 685                         Int_t ichmb = fGeom->GetChamber(xyz0[2], iplane);
 686                         if(ichmb<0) ichmb = fGeom->GetChamber(xyz1[2], iplane);
 687                         if(ichmb<0){
 688                                 // here we should decide what to do with the track. The space between the pads in 2 chambers is 4cm+. Is it making sense to continue building the tracklet here TODO????
 689                                 AliWarning(Form("Track prolongated in the interspace between TRD detectors in plane %d. Skip plane. To be fixed !", iplane));
 690                                 continue;
 691                         }
 692
 693                         // temporary until the functionalities of AliTRDpropagationLayer and AliTRDstackLayer are merged TODO
 694                         AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, ichmb);
 695                         Int_t nrows = pp->GetNrows();
 696                         Double_t stacklength = pp->GetRow0ROC() - pp->GetRowEndROC();/*(nrows - 2) * pp->GetLengthIPad()        + 2 * pp->GetLengthOPad() + (nrows - 1) * pp->GetRowSpacing();*/
 697                         Double_t z0  = fGeom->GetRow0(iplane, ichmb, 0);
 698
 699                         Int_t nClustersChmb = 0;
 700                         AliTRDstackLayer stackLayer[35];
 701                         for(int itb=0; itb<fTimeBinsPerPlane; itb++){
 702                                 const AliTRDpropagationLayer ksmLayer(*(fTrSec[sector]->GetLayer(row1 - itb)));
 703                                 stackLayer[itb] = ksmLayer;
 704 #ifdef DEBUG
 705                                 stackLayer[itb].SetDebugStream(fDebugStreamer);
 706 #endif
 707                                 stackLayer[itb].SetRange(z0 - stacklength, stacklength);
 708                                 stackLayer[itb].SetSector(sector);
 709                                 stackLayer[itb].SetStackNr(ichmb);
 710                                 stackLayer[itb].SetNRows(nrows);
 711                                 stackLayer[itb].SetRecoParam(fRecoParam);
 712                                 stackLayer[itb].BuildIndices();
 713                                 nClustersChmb += stackLayer[itb].GetNClusters();
 714                         }
 715                         if(!nClustersChmb) continue;
 716                         //AliInfo(Form("Detector p[%d] c[%d]. Building tracklet from %d clusters ... ", iplane, ichmb, nClustersChmb));
 717
 718                         tracklet.SetRecoParam(fRecoParam);
 719                         tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
 720                         tracklet.SetPadLength(pp->GetLengthIPad());
 721                         tracklet.SetPlane(iplane);
 722                         //Int_t tbRange   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
 723                         //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
 724                         //tracklet.SetNTimeBinsRange(tbRange);
 725                         tracklet.SetX0(x0);
 726                         tracklet.Init(&t);
 727                         if(!tracklet.AttachClustersIter(stackLayer, 1000.)) continue;
 728                         tracklet.Init(&t);
 729
 730                         //if(!tracklet.AttachClusters(stackLayer, kTRUE)) continue;
 731                         //if(!tracklet.Fit()) continue;
 732                 }
 733                 Int_t ncl = tracklet.GetN();
 734                 //AliInfo(Form("N clusters %d", ncl));
 735
 736                 // Discard tracklet if bad quality.
 737                 //Check if this condition is not already checked during building of the tracklet
 738     if(ncl < fTimeBinsPerPlane * fRecoParam->GetFindableClusters()){
 739                         //AliInfo(Form("Discard tracklet for %d nclusters", ncl));
 740                         continue;
 741                 }
 742
 743                 // load tracklet to the tracker and the track
 744                 Int_t index = SetTracklet(&tracklet);
 745                 t.SetTracklet(&tracklet, iplane, index);
 746
 747                 // Calculate the mean material budget along the path inside the chamber
 748     Double_t param[7];
 749                 AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
 750     // The mean propagation parameters
 751     Double_t xrho = param[0]*param[4]; // density*length
 752     Double_t xx0  = param[1]; // radiation length
 753
 754                 // Propagate and update track
 755                 t.PropagateTo(tracklet.GetX0(), xx0, xrho);
 756           if (!AdjustSector(&t)) break;
 757                 Double_t maxChi2 = t.GetPredictedChi2(&tracklet);
 758                 if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){
 759                         nClustersExpected += ncl;
 760                 }
 761                 // Reset material budget if 2 consecutive gold
 762                 if(iplane>0 && ncl + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
 763
 764                 // Make backup of the track until is gold
 765                 // TO DO update quality check of the track.
 766                 // consider comparison with fTimeBinsRange
 767                 Float_t ratio0 = ncl / Float_t(fTimeBinsPerPlane);
 768                 //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);
 769     //printf("tracklet.GetChi2() %f     [< 18.0]\n", tracklet.GetChi2());
 770                 //printf("ratio0    %f              [>   0.8]\n", ratio0);
 771                 //printf("ratio1     %f             [>   0.6]\n", ratio1);
 772                 //printf("ratio0+ratio1 %f          [>   1.5]\n", ratio0+ratio1);
 773                 //printf("t.GetNCross()  %d         [==    0]\n", t.GetNCross());
 774                 //printf("TMath::Abs(t.GetSnp()) %f [<  0.85]\n", TMath::Abs(t.GetSnp()));
 775                 //printf("t.GetNumberOfClusters() %d [>    20]\n", t.GetNumberOfClusters());
 776
 777                 if (//(tracklet.GetChi2()      <  18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update
 778         (ratio0                  >   0.8) &&
 779         //(ratio1                  >   0.6) &&
 780         //(ratio0+ratio1           >   1.5) &&
 781         (t.GetNCross()           ==    0) &&
 782         (TMath::Abs(t.GetSnp())  <  0.85) &&
 783         (t.GetNumberOfClusters() >    20)) t.MakeBackupTrack();
 784
 785         } // end planes loop
 786
 787 #ifdef DEBUG
 788         if(AliTRDReconstructor::StreamLevel() > 1){
 789                 TTreeSRedirector &cstreamer = *fDebugStreamer;
 790                 cstreamer << "FollowBackProlongation"
 791                         << "ncl="      << nClustersExpected
 792                         << "track.="   << &t
 793                         << "\n";
 794         }
 795 #endif
 796
 797         return nClustersExpected;
 798 }
 799
 800 //____________________________________________________________________
 801 void AliTRDtrackerV1::UnloadClusters()
 802 {
 803   //
 804   // Clears the arrays of clusters and tracks. Resets sectors and timebins
 805   //
 806
 807   Int_t i;
 808   Int_t nentr;
 809
 810   nentr = fClusters->GetEntriesFast();
 811   for (i = 0; i < nentr; i++) {
 812     delete fClusters->RemoveAt(i);
 813   }
 814   fNclusters = 0;
 815
 816         if(fTracklets){
 817         for (i = 0; i < fTracklets->GetEntriesFast(); i++) delete fTracklets->RemoveAt(i);
 818         }
 819
 820   nentr = fSeeds->GetEntriesFast();
 821   for (i = 0; i < nentr; i++) {
 822     delete fSeeds->RemoveAt(i);
 823   }
 824
 825   nentr = fTracks->GetEntriesFast();
 826         for (i = 0; i < nentr; i++) {
 827     delete fTracks->RemoveAt(i);
 828   }
 829
 830   Int_t nsec = AliTRDgeometry::kNsect;
 831   for (i = 0; i < nsec; i++) {
 832     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
 833       fTrSec[i]->GetLayer(pl)->Clear();
 834     }
 835   }
 836
 837 }
 838
 839 //____________________________________________________________________
 840 AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t &idx)
 841 {
 842 // Find tracklet for TRD track <track>
 843 // Parameters
 844 // - track
 845 // - sector
 846 // - plane
 847 // - index
 848 // Output
 849 // tracklet
 850 // index
 851 // Detailed description
 852 //
 853         idx = track->GetTrackletIndex(p);
 854         //AliInfo(Form("looking for tracklet in plane %d idx %d [%d]", p, idx, track->GetTrackletIndex(p)));
 855         AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
 856         //AliInfo(Form("found 0x%x @ %d", tracklet, idx));
 857
 858 //   Int_t *index = track->GetTrackletIndexes();
 859 //   for (UInt_t i = 0; i < 6; i++) AliInfo(Form("index[%d] = %d", i, index[i]));
 860 //
 861 //      for (UInt_t i = 0; i < 6/*kMaxTimeBinIndex*/; i++) {
 862 //     if (index[i] < 0) continue;
 863 //
 864 //              tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index[i]);
 865 //              if(!tracklet) break;
 866 //
 867 //              if(tracklet->GetPlane() != p) continue;
 868 //
 869 //              idx = index[i];
 870 //      }
 871
 872         return tracklet;
 873 }
 874
 875 //____________________________________________________________________
 876 Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet)
 877 {
 878 // Add this tracklet to the list of tracklets stored in the tracker
 879 //
 880 // Parameters
 881 //   - tracklet : pointer to the tracklet to be added to the list
 882 //
 883 // Output
 884 //   - the index of the new tracklet in the tracker tracklets list
 885 //
 886 // Detailed description
 887 // Build the tracklets list if it is not yet created (late initialization)
 888 // and adds the new tracklet to the list.
 889 //
 890         if(!fTracklets){
 891                 fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsect()*kMaxTracksStack);
 892                 fTracklets->SetOwner(kTRUE);
 893         }
 894         Int_t nentries = fTracklets->GetEntriesFast();
 895         new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
 896         return nentries;
 897 }
 898
 899 //____________________________________________________________________
 900 Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
 901                                        , AliESDEvent *esd)
 902 {
 903   //
 904   // Steer tracking for one SM.
 905   //
 906   // Parameters :
 907   //   sector  : Array of (SM) propagation layers containing clusters
 908   //   esd     : The current ESD event. On output it contains the also
 909   //             the ESD (TRD) tracks found in this SM.
 910   //
 911   // Output :
 912   //   Number of tracks found in this TRD supermodule.
 913   //
 914   // Detailed description
 915   //
 916   // 1. Unpack AliTRDpropagationLayers objects for each stack.
 917   // 2. Launch stack tracking.
 918   //    See AliTRDtrackerV1::Clusters2TracksStack() for details.
 919   // 3. Pack results in the ESD event.
 920   //
 921
 922         AliTRDpadPlane *pp = 0x0;
 923
 924         // allocate space for esd tracks in this SM
 925         TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
 926         esdTrackList.SetOwner();
 927         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 928         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 929         const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins);
 930
 931         Int_t ntracks = 0;
 932         Int_t nClStack = 0;
 933         for(int istack = 0; istack<AliTRDpropagationLayer::kZones; istack++){
 934                 AliTRDstackLayer stackLayer[kNPlanes*kNTimeBins];
 935
 936                 nClStack = 0;
 937                 //AliInfo(Form("Processing stack %i ...",istack));
 938                 //AliInfo("Building stack propagation layers ...");
 939                 for(int ilayer=0; ilayer<kNPlanes*nTimeBins; ilayer++){
 940                         pp = fGeom->GetPadPlane((Int_t)(ilayer/nTimeBins), istack);
 941                         Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad()
 942                                              + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
 943                         //Debug
 944                         Double_t z0  = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
 945                         const AliTRDpropagationLayer ksmLayer(*(sector->GetLayer(ilayer)));
 946                         stackLayer[ilayer] = ksmLayer;
 947 #ifdef DEBUG
 948                         stackLayer[ilayer].SetDebugStream(fDebugStreamer);
 949 #endif
 950                         stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
 951                         stackLayer[ilayer].SetSector(sector->GetSector());
 952                         stackLayer[ilayer].SetStackNr(istack);
 953                         stackLayer[ilayer].SetNRows(pp->GetNrows());
 954                         stackLayer[ilayer].SetRecoParam(fRecoParam);
 955                         stackLayer[ilayer].BuildIndices();
 956                         nClStack += stackLayer[ilayer].GetNClusters();
 957                 }
 958                 //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
 959                 if(nClStack < kFindable) continue;
 960                 ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
 961         }
 962         //AliInfo(Form("Found %d tracks in SM", ntracks));
 963
 964         for(int itrack=0; itrack<ntracks; itrack++)
 965           esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
 966
 967         return ntracks;
 968 }
 969
 970 //____________________________________________________________________
 971 Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
 972                                           , TClonesArray *esdTrackList)
 973 {
 974   //
 975   // Make tracks in one TRD stack.
 976   //
 977   // Parameters :
 978   //   layer  : Array of stack propagation layers containing clusters
 979   //   esdTrackList  : Array of ESD tracks found by the stand alone tracker.
 980   //                   On exit the tracks found in this stack are appended.
 981   //
 982   // Output :
 983   //   Number of tracks found in this stack.
 984   //
 985   // Detailed description
 986   //
 987   // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
 988   // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations.
 989   //    See AliTRDtrackerV1::MakeSeeds() for more details.
 990   // 3. Arrange track candidates in decreasing order of their quality
 991   // 4. Classify tracks in 5 categories according to:
 992   //    a) number of layers crossed
 993   //    b) track quality
 994   // 5. Sign clusters by tracks in decreasing order of track quality
 995   // 6. Build AliTRDtrack out of seeding tracklets
 996   // 7. Cook MC label
 997   // 8. Build ESD track and register it to the output list
 998   //
 999
1000         AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
1001         Int_t pars[4]; // MakeSeeds parameters
1002
1003         //Double_t alpha = AliTRDgeometry::GetAlpha();
1004         //Double_t shift = .5 * alpha;
1005         Int_t configs[kNConfigs];
1006
1007         // Build initial seeding configurations
1008         Double_t quality = BuildSeedingConfigs(layer, configs);
1009 #ifdef DEBUG
1010                 if(AliTRDReconstructor::StreamLevel() > 1)
1011                   AliInfo(Form("Plane config %d %d %d Quality %f"
1012                               , configs[0], configs[1], configs[2], quality));
1013 #endif
1014
1015         // Initialize contors
1016         Int_t ntracks,      // number of TRD track candidates
1017               ntracks1,     // number of registered TRD tracks/iter
1018               ntracks2 = 0; // number of all registered TRD tracks in stack
1019         fSieveSeeding = 0;
1020         do{
1021                 // Loop over seeding configurations
1022                 ntracks = 0; ntracks1 = 0;
1023                 for (Int_t iconf = 0; iconf<3; iconf++) {
1024                         pars[0] = configs[iconf];
1025                         pars[1] = layer->GetStackNr();
1026                         pars[2] = ntracks;
1027                         ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars);
1028                         if(ntracks == kMaxTracksStack) break;
1029                 }
1030 #ifdef DEBUG
1031                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in stack %d iteration %d.", ntracks, pars[1], fSieveSeeding));
1032 #endif
1033                 if(!ntracks) break;
1034
1035                 // Sort the seeds according to their quality
1036                 Int_t sort[kMaxTracksStack];
1037                 TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
1038
1039                 // Initialize number of tracks so far and logic switches
1040                 Int_t ntracks0 = esdTrackList->GetEntriesFast();
1041                 Bool_t signedTrack[kMaxTracksStack];
1042                 Bool_t fakeTrack[kMaxTracksStack];
1043                 for (Int_t i=0; i<ntracks; i++){
1044                         signedTrack[i] = kFALSE;
1045                         fakeTrack[i] = kFALSE;
1046                 }
1047                 //AliInfo("Selecting track candidates ...");
1048
1049                 // Sieve clusters in decreasing order of track quality
1050                 Double_t trackParams[7];
1051 //              AliTRDseedV1 *lseed = 0x0;
1052                 Int_t jSieve = 0, candidates;
1053                 do{
1054                         //AliInfo(Form("\t\tITER = %i ", jSieve));
1055
1056                         // Check track candidates
1057                         candidates = 0;
1058                         for (Int_t itrack = 0; itrack < ntracks; itrack++) {
1059                                 Int_t trackIndex = sort[itrack];
1060                                 if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
1061
1062
1063                                 // Calculate track parameters from tracklets seeds
1064                                 Int_t labelsall[1000];
1065                                 Int_t nlabelsall = 0;
1066                                 Int_t naccepted  = 0;
1067                                 Int_t ncl        = 0;
1068                                 Int_t nused      = 0;
1069                                 Int_t nlayers    = 0;
1070                                 Int_t findable   = 0;
1071                                 for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
1072                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
1073                                         if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15)
1074                                           findable++;
1075
1076                                         if(!sseed[jseed].IsOK()) continue;
1077                                         sseed[jseed].UpdateUsed();
1078                                         ncl   += sseed[jseed].GetN2();
1079                                         nused += sseed[jseed].GetNUsed();
1080                                         nlayers++;
1081
1082                                         // Cooking label
1083                                         for (Int_t itime = 0; itime < fTimeBinsPerPlane; itime++) {
1084                                                 if(!sseed[jseed].IsUsable(itime)) continue;
1085                                                 naccepted++;
1086                                                 Int_t tindex = 0, ilab = 0;
1087                                                 while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
1088                                                         labelsall[nlabelsall++] = tindex;
1089                                                         ilab++;
1090                                                 }
1091                                         }
1092                                 }
1093                                 // Filter duplicated tracks
1094                                 if (nused > 30){
1095                                         printf("Skip %d nused %d\n", trackIndex, nused);
1096                                         fakeTrack[trackIndex] = kTRUE;
1097                                         continue;
1098                                 }
1099                                 if (Float_t(nused)/ncl >= .25){
1100                                         printf("Skip %d nused/ncl >= .25\n", trackIndex);
1101                                         fakeTrack[trackIndex] = kTRUE;
1102                                         continue;
1103                                 }
1104
1105                                 // Classify tracks
1106                                 Bool_t skip = kFALSE;
1107                                 switch(jSieve){
1108                                 case 0:
1109                                         if(nlayers < 6) {skip = kTRUE; break;}
1110                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1111                                         break;
1112
1113                                 case 1:
1114                                         if(nlayers < findable){skip = kTRUE; break;}
1115                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
1116                                         break;
1117
1118                                 case 2:
1119                                         if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
1120                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
1121                                         break;
1122
1123                                 case 3:
1124                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
1125                                         break;
1126
1127                                 case 4:
1128                                         if (nlayers == 3){skip = kTRUE; break;}
1129                                         //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
1130                                         break;
1131                                 }
1132                                 if(skip){
1133                                         candidates++;
1134                                         printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
1135                                         continue;
1136                                 }
1137                                 signedTrack[trackIndex] = kTRUE;
1138
1139
1140                                 // Build track label - what happens if measured data ???
1141                                 Int_t labels[1000];
1142                                 Int_t outlab[1000];
1143                                 Int_t nlab = 0;
1144                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1145                                         Int_t jseed = kNPlanes*trackIndex+iLayer;
1146                                         if(!sseed[jseed].IsOK()) continue;
1147                                         for(int ilab=0; ilab<2; ilab++){
1148                                                 if(sseed[jseed].GetLabels(ilab) < 0) continue;
1149                                                 labels[nlab] = sseed[jseed].GetLabels(ilab);
1150                                                 nlab++;
1151                                         }
1152                                 }
1153                                 Freq(nlab,labels,outlab,kFALSE);
1154                                 Int_t   label     = outlab[0];
1155                                 Int_t   frequency = outlab[1];
1156                                 Freq(nlabelsall,labelsall,outlab,kFALSE);
1157                                 Int_t   label1    = outlab[0];
1158                                 Int_t   label2    = outlab[2];
1159                                 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
1160
1161
1162                                 // Sign clusters
1163                                 AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
1164                                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1165                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
1166                                         if(!sseed[jseed].IsOK()) continue;
1167                                         if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
1168                                         sseed[jseed].UseClusters();
1169                                         if(!cl){
1170                                                 Int_t ic = 0;
1171                                                 while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
1172                                                 clusterIndex =  sseed[jseed].GetIndexes(ic);
1173                                         }
1174                                 }
1175                                 if(!cl) continue;
1176
1177
1178                                 // Build track parameters
1179                                 AliTRDseedV1 *lseed =&sseed[trackIndex*6];
1180                                 Int_t idx = 0;
1181                                 while(idx<3 && !lseed->IsOK()) {
1182                                         idx++;
1183                                         lseed++;
1184                                 }
1185                                 Double_t cR = lseed->GetC();
1186                                 trackParams[1] = lseed->GetYref(0);
1187                                 trackParams[2] = lseed->GetZref(0);
1188                                 trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
1189                                 trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
1190                                 trackParams[5] = cR;
1191                                 trackParams[0] = lseed->GetX0();
1192                                 trackParams[6] = layer[0].GetSector();/* *alpha+shift;  // Supermodule*/
1193
1194 #ifdef DEBUG
1195                                 if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]);
1196
1197                                 if(AliTRDReconstructor::StreamLevel() >= 1){
1198                                         Int_t sector = layer[0].GetSector();
1199                                         Int_t nclusters = 0;
1200                                         AliTRDseedV1 *dseed[6];
1201                                         for(int is=0; is<6; is++){
1202                                                 dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is]);
1203                                                 dseed[is]->SetOwner();
1204                                                 nclusters += sseed[is].GetN2();
1205                                                 //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
1206                                         }
1207                                         //Int_t eventNrInFile = esd->GetEventNumberInFile();
1208                                         //AliInfo(Form("Number of clusters %d.", nclusters));
1209                                         TTreeSRedirector &cstreamer = *fDebugStreamer;
1210                                         cstreamer << "Clusters2TracksStack"
1211                                                 << "Iter="      << fSieveSeeding
1212                                                 << "Like="      << fTrackQuality[trackIndex]
1213                                                 << "S0.="       << dseed[0]
1214                                                 << "S1.="       << dseed[1]
1215                                                 << "S2.="       << dseed[2]
1216                                                 << "S3.="       << dseed[3]
1217                                                 << "S4.="       << dseed[4]
1218                                                 << "S5.="       << dseed[5]
1219                                                 << "p0=" << trackParams[0]
1220                                                 << "p1=" << trackParams[1]
1221                                                 << "p2=" << trackParams[2]
1222                                                 << "p3=" << trackParams[3]
1223                                                 << "p4=" << trackParams[4]
1224                                                 << "p5=" << trackParams[5]
1225                                                 << "p6=" << trackParams[6]
1226                                                 << "Sector="    << sector
1227                                                 << "Stack="     << pars[1]
1228                                                 << "Label="     << label
1229                                                 << "Label1="    << label1
1230                                                 << "Label2="    << label2
1231                                                 << "FakeRatio=" << fakeratio
1232                                                 << "Freq="      << frequency
1233                                                 << "Ncl="       << ncl
1234                                                 << "NLayers="   << nlayers
1235                                                 << "Findable="  << findable
1236                                                 << "NUsed="     << nused
1237                                                 << "\n";
1238                                         //???for(int is=0; is<6; is++) delete dseed[is];
1239                                 }
1240 #endif
1241
1242                                 AliTRDtrackV1 *track = AliTRDtrackerV1::MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
1243                                 if(!track){
1244                                         AliWarning("Fail to build a TRD Track.");
1245                                         continue;
1246                                 }
1247                                 AliInfo("End of MakeTrack()");
1248                                 AliESDtrack esdTrack;
1249                                 esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
1250                                 esdTrack.SetLabel(track->GetLabel());
1251                                 new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
1252                                 ntracks1++;
1253                         }
1254
1255                         jSieve++;
1256                 } while(jSieve<5 && candidates); // end track candidates sieve
1257                 if(!ntracks1) break;
1258
1259                 // increment counters
1260                 ntracks2 += ntracks1;
1261                 fSieveSeeding++;
1262
1263                 // Rebuild plane configurations and indices taking only unused clusters into account
1264                 quality = BuildSeedingConfigs(layer, configs);
1265                 //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
1266
1267                 for(Int_t il = 0; il < kNPlanes * fTimeBinsPerPlane; il++) layer[il].BuildIndices(fSieveSeeding);
1268
1269 #ifdef DEBUG
1270                                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
1271 #endif
1272         } while(fSieveSeeding<10); // end stack clusters sieve
1273
1274
1275
1276         //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
1277
1278         return ntracks2;
1279 }
1280
1281 //___________________________________________________________________
1282 Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
1283                                             , Int_t *configs)
1284 {
1285   //
1286   // Assign probabilities to chambers according to their
1287   // capability of producing seeds.
1288   //
1289   // Parameters :
1290   //
1291   //   layers : Array of stack propagation layers for all 6 chambers in one stack
1292   //   configs : On exit array of configuration indexes (see GetSeedingConfig()
1293   // for details) in the decreasing order of their seeding probabilities.
1294   //
1295   // Output :
1296   //
1297   //  Return top configuration quality
1298   //
1299   // Detailed description:
1300   //
1301   // To each chamber seeding configuration (see GetSeedingConfig() for
1302   // the list of all configurations) one defines 2 quality factors:
1303   //  - an apriori topological quality (see GetSeedingConfig() for details) and
1304   //  - a data quality based on the uniformity of the distribution of
1305   //    clusters over the x range (time bins population). See CookChamberQA() for details.
1306   // The overall chamber quality is given by the product of this 2 contributions.
1307   //
1308
1309         Double_t chamberQA[kNPlanes];
1310         for(int iplane=0; iplane<kNPlanes; iplane++){
1311                 chamberQA[iplane] = MakeSeedingPlanes(&layers[iplane*fTimeBinsPerPlane]);
1312                 //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
1313         }
1314
1315         Double_t tconfig[kNConfigs];
1316         Int_t planes[4];
1317         for(int iconf=0; iconf<kNConfigs; iconf++){
1318                 GetSeedingConfig(iconf, planes);
1319                 tconfig[iconf] = fgTopologicQA[iconf];
1320                 for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]];
1321         }
1322
1323         TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
1324         return tconfig[configs[0]];
1325 }
1326
1327 //____________________________________________________________________
1328 Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
1329                                , AliTRDseedV1 *sseed
1330                                , Int_t *ipar)
1331 {
1332   //
1333   // Make tracklet seeds in the TRD stack.
1334   //
1335   // Parameters :
1336   //   layers : Array of stack propagation layers containing clusters
1337   //   sseed  : Array of empty tracklet seeds. On exit they are filled.
1338   //   ipar   : Control parameters:
1339   //       ipar[0] -> seeding chambers configuration
1340   //       ipar[1] -> stack index
1341   //       ipar[2] -> number of track candidates found so far
1342   //
1343   // Output :
1344   //   Number of tracks candidates found.
1345   //
1346   // Detailed description
1347   //
1348   // The following steps are performed:
1349   // 1. Select seeding layers from seeding chambers
1350   // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
1351   //   The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
1352   //   this order. The parameters controling the range of accepted clusters in
1353   //   layer 0, 1, and 2 are defined in AliTRDstackLayer::BuildCond().
1354   // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
1355   // 4. Initialize seeding tracklets in the seeding chambers.
1356   // 5. Filter 0.
1357   //   Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
1358   //   Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
1359   // 6. Attach clusters to seeding tracklets and find linear approximation of
1360   //   the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
1361   //   clusters used by current seeds should not exceed ... (25).
1362   // 7. Filter 1.
1363   //   All 4 seeding tracklets should be correctly constructed (see
1364   //   AliTRDseedV1::AttachClustersIter())
1365   // 8. Helix fit of the seeding tracklets
1366   // 9. Filter 2.
1367   //   Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
1368   // 10. Extrapolation of the helix fit to the other 2 chambers:
1369   //    a) Initialization of extrapolation tracklet with fit parameters
1370   //    b) Helix fit of tracklets
1371   //    c) Attach clusters and linear interpolation to extrapolated tracklets
1372   //    d) Helix fit of tracklets
1373   // 11. Improve seeding tracklets quality by reassigning clusters.
1374   //      See AliTRDtrackerV1::ImproveSeedQuality() for details.
1375   // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
1376   // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
1377   // 14. Cooking labels for tracklets. Should be done only for MC
1378   // 15. Register seeds.
1379   //
1380
1381         AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
1382         AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
1383         Int_t ncl, mcl; // working variable for looping over clusters
1384         Int_t index[AliTRDstackLayer::kMaxClustersLayer], jndex[AliTRDstackLayer::kMaxClustersLayer];
1385         // chi2 storage
1386         // chi2[0] = tracklet chi2 on the Z direction
1387         // chi2[1] = tracklet chi2 on the R direction
1388         Double_t chi2[4];
1389
1390
1391         // this should be data member of AliTRDtrack
1392         Double_t seedQuality[kMaxTracksStack];
1393
1394         // unpack control parameters
1395         Int_t config  = ipar[0];
1396         Int_t istack  = ipar[1];
1397         Int_t ntracks = ipar[2];
1398         Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);
1399 #ifdef DEBUG
1400                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
1401 #endif
1402
1403         // Init chambers geometry
1404         Int_t det/*, tbRange[6]*/; // time bins inside the detector geometry
1405         Double_t hL[kNPlanes];       // Tilting angle
1406         Float_t padlength[kNPlanes]; // pad lenghts
1407         AliTRDpadPlane *pp;
1408         for(int il=0; il<kNPlanes; il++){
1409                 pp = fGeom->GetPadPlane(il, istack); // istack has to be imported
1410                 hL[il]        = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
1411                 padlength[il] = pp->GetLengthIPad();
1412                 det           = il; // to be fixed !!!!!
1413                 //tbRange[il]   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
1414                 //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
1415         }
1416
1417         Double_t cond0[4], cond1[4], cond2[4];
1418         // make seeding layers (to be moved in Clusters2TracksStack)
1419         AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0};
1420         for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * fTimeBinsPerPlane], planes[isl]);
1421
1422
1423         // Start finding seeds
1424         Int_t icl = 0;
1425         while((c[3] = (*layer[3])[icl++])){
1426                 if(!c[3]) continue;
1427                 layer[0]->BuildCond(c[3], cond0, 0);
1428                 layer[0]->GetClusters(cond0, index, ncl);
1429                 Int_t jcl = 0;
1430                 while(jcl<ncl) {
1431                         c[0] = (*layer[0])[index[jcl++]];
1432                         if(!c[0]) continue;
1433                         Double_t dx    = c[3]->GetX() - c[0]->GetX();
1434                         Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
1435                         Double_t phi   = (c[3]->GetY() - c[0]->GetY())/dx;
1436                         layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
1437                         layer[1]->GetClusters(cond1, jndex, mcl);
1438
1439                         Int_t kcl = 0;
1440                         while(kcl<mcl) {
1441                                 c[1] = (*layer[1])[jndex[kcl++]];
1442                                 if(!c[1]) continue;
1443                                 layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
1444                                 c[2] = layer[2]->GetNearestCluster(cond2);
1445                                 if(!c[2]) continue;
1446
1447                                 //AliInfo("Seeding clusters found. Building seeds ...");
1448                                 //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
1449                                 for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
1450
1451                                 fFitter->Reset();
1452
1453                                 fFitter->FitRieman(c, kNSeedPlanes);
1454
1455                                 chi2[0] = 0.; chi2[1] = 0.;
1456                                 AliTRDseedV1 *tseed = 0x0;
1457                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
1458                                         Int_t jLayer = planes[iLayer];
1459                                         tseed = &cseed[jLayer];
1460                                         tseed->SetRecoParam(fRecoParam);
1461                                         tseed->SetPlane(jLayer);
1462                                         tseed->SetTilt(hL[jLayer]);
1463                                         tseed->SetPadLength(padlength[jLayer]);
1464                                         //tseed->SetNTimeBinsRange(tbRange[jLayer]);
1465                                         tseed->SetX0(layer[iLayer]->GetX());//layers[jLayer*fTimeBinsPerPlane].GetX());
1466
1467                                         tseed->Init(fFitter->GetRiemanFitter());
1468                                         // temporary until new AttachClusters()
1469                                         tseed->SetX0(layers[(jLayer+1)*fTimeBinsPerPlane-1].GetX());
1470                                         chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
1471                                         chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
1472                                 }
1473
1474                                 Bool_t isFake = kFALSE;
1475                                 if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1476                                 if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1477                                 if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
1478 #ifdef DEBUG
1479                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1480                                         Float_t yref[4], ycluster[4];
1481                                         for(int il=0; il<4; il++){
1482                                                 tseed = &cseed[planes[il]];
1483                                                 yref[il] = tseed->GetYref(0);
1484                                                 ycluster[il] = c[il]->GetY();
1485                                         }
1486                                         Float_t threshold = .5;//1./(3. - sLayer);
1487                                         Int_t ll = c[3]->GetLabel(0);
1488                                         TTreeSRedirector &cs0 = *fDebugStreamer;
1489                                                         cs0 << "MakeSeeds0"
1490                                                         <<"isFake=" << isFake
1491                                                         <<"label=" << ll
1492                                                         <<"threshold=" << threshold
1493                                                         <<"chi2=" << chi2[1]
1494                                                         <<"yref0="<<yref[0]
1495                                                         <<"yref1="<<yref[1]
1496                                                         <<"yref2="<<yref[2]
1497                                                         <<"yref3="<<yref[3]
1498                                                         <<"ycluster0="<<ycluster[0]
1499                                                         <<"ycluster1="<<ycluster[1]
1500                                                         <<"ycluster2="<<ycluster[2]
1501                                                         <<"ycluster3="<<ycluster[3]
1502                                                         <<"\n";
1503                                 }
1504 #endif
1505
1506                                 if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
1507                                         //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
1508                                         continue;
1509                                 }
1510                                 if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
1511                                         //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
1512                                         continue;
1513                                 }
1514                                 //AliInfo("Passed chi2 filter.");
1515
1516 #ifdef DEBUG
1517                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1518                                         Float_t minmax[2] = { -100.0,  100.0 };
1519                                         for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
1520                                                 Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
1521                                                 if (max < minmax[1]) minmax[1] = max;
1522                                                 Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
1523                                                 if (min > minmax[0]) minmax[0] = min;
1524                                         }
1525                                         Double_t xpos[4];
1526                                         for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
1527                                         TTreeSRedirector &cstreamer = *fDebugStreamer;
1528                                                         cstreamer << "MakeSeeds1"
1529                                                 << "isFake=" << isFake
1530                                                 << "config="   << config
1531                                                 << "Cl0.="   << c[0]
1532                                                 << "Cl1.="   << c[1]
1533                                                 << "Cl2.="   << c[2]
1534                                                 << "Cl3.="   << c[3]
1535                                                 << "X0="     << xpos[0] //layer[sLayer]->GetX()
1536                                                 << "X1="     << xpos[1] //layer[sLayer + 1]->GetX()
1537                                                 << "X2="     << xpos[2] //layer[sLayer + 2]->GetX()
1538                                                 << "X3="     << xpos[3] //layer[sLayer + 3]->GetX()
1539                                                 << "Y2exp="  << cond2[0]
1540                                                 << "Z2exp="  << cond2[1]
1541                                                 << "Chi2R="  << chi2[0]
1542                                                 << "Chi2Z="  << chi2[1]
1543                                                 << "Seed0.=" << &cseed[planes[0]]
1544                                                 << "Seed1.=" << &cseed[planes[1]]
1545                                                 << "Seed2.=" << &cseed[planes[2]]
1546                                                 << "Seed3.=" << &cseed[planes[3]]
1547                                                 << "Zmin="   << minmax[0]
1548                                                 << "Zmax="   << minmax[1]
1549                                                 << "\n" ;
1550                                 }
1551 #endif
1552                                 // try attaching clusters to tracklets
1553                                 Int_t nUsedCl = 0;
1554                                 Int_t nlayers = 0;
1555                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
1556                                         Int_t jLayer = planes[iLayer];
1557                                         if(!cseed[jLayer].AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 5., kFALSE, c[iLayer])) continue;
1558                                         nUsedCl += cseed[jLayer].GetNUsed();
1559                                         if(nUsedCl > 25) break;
1560                                         nlayers++;
1561                                 }
1562                                 if(nlayers < kNSeedPlanes){
1563                                         //AliInfo("Failed updating all seeds.");
1564                                         continue;
1565                                 }
1566                                 // fit tracklets and cook likelihood
1567                                 chi2[0] = 0.; chi2[1] = 0.;
1568                                 fFitter->FitRieman(&cseed[0], &planes[0]);
1569                                 AliRieman *rim = fFitter->GetRiemanFitter();
1570                                 for(int iLayer=0; iLayer<4; iLayer++){
1571                                         cseed[planes[iLayer]].Init(rim);
1572                                         chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
1573                                         chi2[1] += cseed[planes[iLayer]].GetChi2Y();
1574                                 }
1575                                 Double_t chi2r = chi2[1], chi2z = chi2[0];
1576                                 Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
1577                                 if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
1578                                         //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
1579                                         continue;
1580                                 }
1581                                 //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
1582
1583
1584                                 // book preliminary results
1585                                 seedQuality[ntracks] = like;
1586                                 fSeedLayer[ntracks]  = config;/*sLayer;*/
1587
1588                                 // attach clusters to the extrapolation seeds
1589                                 Int_t lextrap[2];
1590                                 GetExtrapolationConfig(config, lextrap);
1591                                 Int_t nusedf   = 0; // debug value
1592                                 for(int iLayer=0; iLayer<2; iLayer++){
1593                                         Int_t jLayer = lextrap[iLayer];
1594
1595                                         // prepare extrapolated seed
1596                                         cseed[jLayer].Reset();
1597                                         cseed[jLayer].SetRecoParam(fRecoParam);
1598                                         cseed[jLayer].SetPlane(jLayer);
1599                                         cseed[jLayer].SetTilt(hL[jLayer]);
1600                                         cseed[jLayer].SetX0(layers[(jLayer +1) * fTimeBinsPerPlane-1].GetX());
1601                                         cseed[jLayer].SetPadLength(padlength[jLayer]);
1602                                         //cseed[jLayer].SetNTimeBinsRange(tbRange[jLayer]);
1603                                         cseed[jLayer].Init(rim);
1604 //                                      AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*fTimeBinsPerPlane], &cseed[jLayer]);
1605 //                                      if(cd == 0x0) continue;
1606
1607                                         // fit extrapolated seed
1608                                         AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
1609                                         if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
1610                                         if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
1611                                         AliTRDseedV1 tseed = cseed[jLayer];
1612                                         if(!tseed.AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 1000.)) continue;
1613                                         cseed[jLayer] = tseed;
1614                                         nusedf += cseed[jLayer].GetNUsed(); // debug value
1615                                         AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
1616                                 }
1617                                 //AliInfo("Extrapolation done.");
1618
1619                                 ImproveSeedQuality(layers, cseed);
1620                                 //AliInfo("Improve seed quality done.");
1621
1622                                 nlayers   = 0;
1623                                 Int_t nclusters = 0;
1624                                 Int_t findable  = 0;
1625                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1626                                         if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++;
1627                                         if (!cseed[iLayer].IsOK()) continue;
1628                                         nclusters += cseed[iLayer].GetN2();
1629                                         nlayers++;
1630                                 }
1631                                 if (nlayers < 3){
1632                                         //AliInfo("Failed quality check on seeds.");
1633                                         continue;
1634                                 }
1635
1636                                 // fit full track and cook likelihoods
1637                                 fFitter->FitRieman(&cseed[0]);
1638                                 Double_t chi2ZF = 0., chi2RF = 0.;
1639                                 for(int ilayer=0; ilayer<6; ilayer++){
1640                                         cseed[ilayer].Init(fFitter->GetRiemanFitter());
1641                                         if (!cseed[ilayer].IsOK()) continue;
1642                                         //tchi2 = cseed[ilayer].GetChi2Z();
1643                                         //printf("layer %d chi2 %e\n", ilayer, tchi2);
1644                                         chi2ZF += cseed[ilayer].GetChi2Z();
1645                                         chi2RF += cseed[ilayer].GetChi2Y();
1646                                 }
1647                                 chi2ZF /= TMath::Max((nlayers - 3.), 1.);
1648                                 chi2RF /= TMath::Max((nlayers - 3.), 1.);
1649
1650                                 // do the final track fitting
1651                                 fFitter->SetLayers(nlayers);
1652 #ifdef DEBUG
1653                                 fFitter->SetDebugStream(fDebugStreamer);
1654 #endif
1655                                 fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
1656                                 Double_t param[3];
1657                                 Double_t chi2[2];
1658                                 fFitter->GetHyperplaneFitResults(param);
1659                                 fFitter->GetHyperplaneFitChi2(chi2);
1660                                 //AliInfo("Hyperplane fit done\n");
1661
1662                                 // finalize tracklets
1663                                 Int_t labels[12];
1664                                 Int_t outlab[24];
1665                                 Int_t nlab = 0;
1666                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1667                                         if (!cseed[iLayer].IsOK()) continue;
1668
1669                                         if (cseed[iLayer].GetLabels(0) >= 0) {
1670                                                 labels[nlab] = cseed[iLayer].GetLabels(0);
1671                                                 nlab++;
1672                                         }
1673
1674                                         if (cseed[iLayer].GetLabels(1) >= 0) {
1675                                                 labels[nlab] = cseed[iLayer].GetLabels(1);
1676                                                 nlab++;
1677                                         }
1678                                 }
1679                                 Freq(nlab,labels,outlab,kFALSE);
1680                                 Int_t label     = outlab[0];
1681                                 Int_t frequency = outlab[1];
1682                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
1683                                         cseed[iLayer].SetFreq(frequency);
1684                                         cseed[iLayer].SetC(param[1]/*cR*/);
1685                                         cseed[iLayer].SetCC(param[0]/*cC*/);
1686                                         cseed[iLayer].SetChi2(chi2[0]);
1687                                         cseed[iLayer].SetChi2Z(chi2ZF);
1688                                 }
1689
1690 #ifdef DEBUG
1691                                 if(AliTRDReconstructor::StreamLevel() >= 2){
1692                                         Double_t curv = (fFitter->GetRiemanFitter())->GetC();
1693                                         TTreeSRedirector &cstreamer = *fDebugStreamer;
1694                                         cstreamer << "MakeSeeds2"
1695                                                 << "C="       << curv
1696                                                 << "Chi2R="   << chi2r
1697                                                 << "Chi2Z="   << chi2z
1698                                                 << "Chi2TR="  << chi2[0]
1699                                                 << "Chi2TC="  << chi2[1]
1700                                                 << "Chi2RF="  << chi2RF
1701                                                 << "Chi2ZF="  << chi2ZF
1702                                                 << "Ncl="     << nclusters
1703                                                 << "Nlayers=" << nlayers
1704                                                 << "NUsedS="  << nUsedCl
1705                                                 << "NUsed="   << nusedf
1706                                                 << "Findable" << findable
1707                                                 << "Like="    << like
1708                                                 << "S0.="     << &cseed[0]
1709                                                 << "S1.="     << &cseed[1]
1710                                                 << "S2.="     << &cseed[2]
1711                                                 << "S3.="     << &cseed[3]
1712                                                 << "S4.="     << &cseed[4]
1713                                                 << "S5.="     << &cseed[5]
1714                                                 << "Label="   << label
1715                                                 << "Freq="    << frequency
1716                                                 << "\n";
1717                                 }
1718 #endif
1719
1720                                 ntracks++;
1721                                 if(ntracks == kMaxTracksStack){
1722                                         AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
1723                                         return ntracks;
1724                                 }
1725                                 cseed += 6;
1726                         }
1727                 }
1728         }
1729         for(int isl=0; isl<4; isl++) delete layer[isl];
1730
1731         return ntracks;
1732 }
1733
1734 //_____________________________________________________________________________
1735 AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
1736 {
1737   //
1738   // Build a TRD track out of tracklet candidates
1739   //
1740   // Parameters :
1741   //   seeds  : array of tracklets
1742   //   params : track parameters (see MakeSeeds() function body for a detailed description)
1743   //
1744   // Output :
1745   //   The TRD track.
1746   //
1747   // Detailed description
1748   //
1749   // To be discussed with Marian !!
1750   //
1751
1752   Double_t alpha = AliTRDgeometry::GetAlpha();
1753   Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
1754   Double_t c[15];
1755
1756   c[ 0] = 0.2;
1757   c[ 1] = 0.0; c[ 2] = 2.0;
1758   c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
1759   c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0;  c[ 9] = 0.1;
1760   c[10] = 0.0; c[11] = 0.0; c[12] = 0.0;  c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
1761
1762   AliTRDtrackV1 *track = new AliTRDtrackV1(seeds, &params[1], c, params[0], params[6]*alpha+shift);
1763         track->PropagateTo(params[0]-5.0);
1764   track->ResetCovariance(1);
1765   Int_t nc = FollowBackProlongation(*track);
1766         AliInfo(Form("N clusters for track %d", nc));
1767         if (nc < 30) {
1768     delete track;
1769     track = 0x0;
1770   } else {
1771 //     track->CookdEdx();
1772 //     track->CookdEdxTimBin(-1);
1773 //     CookLabel(track, 0.9);
1774   }
1775
1776   return track;
1777 }
1778
1779 //____________________________________________________________________
1780 void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const
1781 {
1782         // to be implemented, preferably at the level of TRD tracklet. !!!!!!!
1783 }
1784
1785 //____________________________________________________________________
1786 void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
1787                                        , AliTRDseedV1 *cseed)
1788 {
1789   //
1790   // Sort tracklets according to "quality" and try to "improve" the first 4 worst
1791   //
1792   // Parameters :
1793   //  layers : Array of propagation layers for a stack/supermodule
1794   //  cseed  : Array of 6 seeding tracklets which has to be improved
1795   //
1796   // Output :
1797   //   cssed : Improved seeds
1798   //
1799   // Detailed description
1800   //
1801   // Iterative procedure in which new clusters are searched for each
1802   // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
1803   // can be maximized. If some optimization is found the old seeds are replaced.
1804   //
1805
1806         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1807         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1808
1809         // make a local working copy
1810         AliTRDseedV1 bseed[6];
1811         for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
1812
1813
1814         Float_t lastquality = 10000.0;
1815         Float_t lastchi2    = 10000.0;
1816         Float_t chi2        =  1000.0;
1817
1818         for (Int_t iter = 0; iter < 4; iter++) {
1819                 Float_t sumquality = 0.0;
1820                 Float_t squality[6];
1821                 Int_t   sortindexes[6];
1822
1823                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1824                         squality[jLayer]  = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
1825                         sumquality +=squality[jLayer];
1826                 }
1827                 if ((sumquality >= lastquality) || (chi2       >     lastchi2)) break;
1828
1829
1830                 lastquality = sumquality;
1831                 lastchi2    = chi2;
1832                 if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
1833
1834
1835                 TMath::Sort(6, squality, sortindexes, kFALSE);
1836                 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1837                         Int_t bLayer = sortindexes[jLayer];
1838                         bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE);
1839                 }
1840
1841                 chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
1842         } // Loop: iter
1843 }
1844
1845 //____________________________________________________________________
1846 Double_t  AliTRDtrackerV1::MakeSeedingPlanes(AliTRDstackLayer *layers)
1847 {
1848   //
1849   // Calculate plane quality for seeding.
1850   //
1851   //
1852   // Parameters :
1853   //   layers : Array of propagation layers for this plane.
1854   //
1855   // Output :
1856   //   plane quality factor for seeding
1857   //
1858   // Detailed description
1859   //
1860   // The quality of the plane for seeding is higher if:
1861   //  1. the average timebin population is closer to an integer number
1862   //  2. the distribution of clusters/timebin is closer to a uniform distribution.
1863   //    - the slope of the first derivative of a parabolic fit is small or
1864   //    - the slope of a linear fit is small
1865   //
1866
1867         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1868         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1869
1870 //      Double_t x;
1871 //      TLinearFitter fitter(1, "pol1");
1872 //      fitter.ClearPoints();
1873         Int_t ncl = 0;
1874         Int_t nused = 0;
1875         Int_t nClLayer;
1876         for(int itb=0; itb<nTimeBins; itb++){
1877                 //x = layer[itb].GetX();
1878                 //printf("x[%d] = %f nCls %d\n", itb, x, layer[itb].GetNClusters());
1879                 //if(!layer[itb].GetNClusters()) continue;
1880                 //fitter.AddPoint(&x, layer[itb].GetNClusters(), 1.);
1881                 nClLayer = layers[itb].GetNClusters();
1882                 ncl += nClLayer;
1883                 for(Int_t incl = 0; incl < nClLayer; incl++)
1884                         if((layers[itb].GetCluster(incl))->IsUsed()) nused++;
1885         }
1886
1887         // calculate the deviation of the mean number of clusters from the
1888         // closest integer values
1889         Float_t nclMed = float(ncl-nused)/nTimeBins;
1890         Int_t ncli = Int_t(nclMed);
1891         Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
1892         nclDev -= (nclDev>.5) && ncli ? 0. : 1.;
1893         /*Double_t quality = */ return TMath::Exp(2.*nclDev);
1894
1895 //      // get slope of the derivative
1896 //      if(!fitter.Eval()) return quality;
1897 //      fitter.PrintResults(3);
1898 //      Double_t a = fitter.GetParameter(1);
1899 //
1900 //      printf("ncl_dev(%f)  a(%f)\n", ncl_dev, a);
1901 //      return quality*TMath::Exp(-a);
1902 }
1903
1904 //____________________________________________________________________
1905 Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
1906                                        , Int_t planes[4]
1907                                        , Double_t *chi2)
1908 {
1909   //
1910   // Calculate the probability of this track candidate.
1911   //
1912   // Parameters :
1913   //   cseeds : array of candidate tracklets
1914   //   planes : array of seeding planes (see seeding configuration)
1915   //   chi2   : chi2 values (on the Z and Y direction) from the rieman fit of the track.
1916   //
1917   // Output :
1918   //   likelihood value
1919   //
1920   // Detailed description
1921   //
1922   // The track quality is estimated based on the following 4 criteria:
1923   //  1. precision of the rieman fit on the Y direction (likea)
1924   //  2. chi2 on the Y direction (likechi2y)
1925   //  3. chi2 on the Z direction (likechi2z)
1926   //  4. number of attached clusters compared to a reference value
1927   //     (see AliTRDrecoParam::fkFindable) (likeN)
1928   //
1929   // The distributions for each type of probabilities are given below as of
1930   // (date). They have to be checked to assure consistency of estimation.
1931   //
1932
1933         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1934         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1935         // ratio of the total number of clusters/track which are expected to be found by the tracker.
1936         Float_t fgFindable = fRecoParam->GetFindableClusters();
1937
1938
1939         Int_t nclusters = 0;
1940         Double_t sumda = 0.;
1941         for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
1942                 Int_t jlayer = planes[ilayer];
1943                 nclusters += cseed[jlayer].GetN2();
1944                 sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
1945         }
1946         Double_t likea     = TMath::Exp(-sumda*10.6);
1947         Double_t likechi2y  = 0.0000000001;
1948         if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
1949         Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
1950         Int_t enc = Int_t(fgFindable*4.*nTimeBins);     // Expected Number Of Clusters, normally 72
1951         Double_t likeN     = TMath::Exp(-(enc - nclusters) * 0.19);
1952
1953         Double_t like      = likea * likechi2y * likechi2z * likeN;
1954
1955 #ifdef DEBUG
1956         //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));
1957         if(AliTRDReconstructor::StreamLevel() >= 2){
1958                 TTreeSRedirector &cstreamer = *fDebugStreamer;
1959                 cstreamer << "CookLikelihood"
1960                         << "sumda="     << sumda
1961                         << "chi0="      << chi2[0]
1962                         << "chi1="      << chi2[1]
1963                         << "likea="     << likea
1964                         << "likechi2y=" << likechi2y
1965                         << "likechi2z=" << likechi2z
1966                         << "nclusters=" << nclusters
1967                         << "likeN="     << likeN
1968                         << "like="      << like
1969                         << "\n";
1970         }
1971 #endif
1972
1973         return like;
1974 }
1975
1976 //___________________________________________________________________
1977 void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
1978                                    , Int_t *planes
1979                                    , Double_t *params)
1980 {
1981   //
1982   // Determines the Mean number of clusters per layer.
1983   // Needed to determine good Seeding Layers
1984   //
1985   // Parameters:
1986   //    - Array of AliTRDstackLayers
1987   //    - Container for the params
1988   //
1989   // Detailed description
1990   //
1991   // Two Iterations:
1992   // In the first Iteration the mean is calculted using all layers.
1993   // After this, all layers outside the 1-sigma-region are rejected.
1994   // Then the mean value and the standard-deviation are calculted a second
1995   // time in order to select all layers in the 1-sigma-region as good-candidates.
1996   //
1997
1998         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1999         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2000
2001         Float_t mean = 0, stdev = 0;
2002         Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
2003         Int_t position = 0;
2004         memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
2005         memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
2006         Int_t nused = 0;
2007         for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
2008                 for(Int_t ils = 0; ils < nTimeBins; ils++){
2009                         position = planes[ipl]*nTimeBins + ils;
2010                         ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters();
2011                         nused = 0;
2012                         for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
2013                                 if((layers[position].GetCluster(icl))->IsUsed()) nused++;
2014                         ncl[ipl * nTimeBins + ils] -= nused;
2015                 }
2016         }
2017         // Declaration of quartils:
2018         //Double_t qvals[3] = {0.0, 0.0, 0.0};
2019         //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
2020         // Iterations
2021         Int_t counter;
2022         Double_t *array;
2023         Int_t *limit;
2024         Int_t nLayers = nTimeBins * kNSeedPlanes;
2025         for(Int_t iter = 0; iter < 2; iter++){
2026                 array = (iter == 0) ? &ncl[0] : &mcl[0];
2027                 limit = (iter == 0) ? &nLayers : &counter;
2028                 counter = 0;
2029                 if(iter == 1){
2030                         for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){
2031                                 if((ncl[i] >  mean + stdev) || (ncl[i] <  mean - stdev)) continue; // Outside 1-sigma region
2032 //                              if((ncl[i] >  qvals[2]) || (ncl[i] <  qvals[0])) continue; // Outside 1-sigma region
2033                                 if(ncl[i] == 0) continue;                                                // 0-Layers also rejected
2034                                 mcl[counter] = ncl[i];
2035                                 counter++;
2036                         }
2037                 }
2038                 if(*limit == 0) break;
2039                 printf("Limit = %d\n", *limit);
2040                 //using quartils instead of mean and RMS
2041 //              TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
2042                 mean = TMath::Median(*limit, array, 0x0);
2043                 stdev  = TMath::RMS(*limit, array);
2044         }
2045 //      printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
2046 //      memcpy(params,qvals,sizeof(Double_t)*3);
2047         params[1] = (Double_t)TMath::Nint(mean);
2048         params[0] = (Double_t)TMath::Nint(mean - stdev);
2049         params[2] = (Double_t)TMath::Nint(mean + stdev);
2050
2051 }
2052
2053 //___________________________________________________________________
2054 Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
2055                                       , Double_t *params)
2056 {
2057   //
2058   // Algorithm to find optimal seeding layer
2059   // Layers inside one sigma region (given by Quantiles) are sorted
2060   // according to their difference.
2061   // All layers outside are sorted according t
2062   //
2063   // Parameters:
2064   //     - Array of AliTRDstackLayers (in the current plane !!!)
2065   //     - Container for the Indices of the seeding Layer candidates
2066   //
2067   // Output:
2068   //     - Number of Layers inside the 1-sigma-region
2069   //
2070   // The optimal seeding layer should contain the mean number of
2071   // custers in the layers in one chamber.
2072   //
2073
2074         //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
2075         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2076         const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer;
2077         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2078         Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
2079         memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
2080         memset(indices, 0, sizeof(Int_t)*kNTimeBins);
2081         memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
2082         Int_t nused = 0;
2083         for(Int_t ils = 0; ils < nTimeBins; ils++){
2084                 ncl[ils] = layers[ils].GetNClusters();
2085                 nused = 0;
2086                 for(Int_t icl = 0; icl < ncl[ils]; icl++)
2087                         if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
2088                 ncl[ils] -= nused;
2089         }
2090
2091         Float_t mean = params[1];
2092         for(Int_t ils = 0; ils < nTimeBins; ils++){
2093                 memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils));
2094                 indices[bins[ncl[ils]+1]] = ils;
2095                 for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
2096                         bins[i]++;
2097         }
2098
2099         //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
2100         Int_t sbin = -1;
2101         Int_t nElements;
2102         Int_t position = 0;
2103         TRandom *r = new TRandom();
2104         Int_t iter = 0;
2105         while(1){
2106                 while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
2107                         // Randomly selecting one bin
2108                         sbin = (Int_t)r->Poisson(mean);
2109                 }
2110                 printf("Bin = %d\n",sbin);
2111                 //Randomly selecting one Layer in the bin
2112                 nElements = bins[sbin + 1] - bins[sbin];
2113                 printf("nElements = %d\n", nElements);
2114                 if(iter == 5){
2115                         position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
2116                         break;
2117                 }
2118                 else if(nElements==0){
2119                         iter++;
2120                         continue;
2121                 }
2122                 position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
2123                 break;
2124         }
2125         delete r;
2126         return indices[position];
2127 }
2128
2129 //____________________________________________________________________
2130 AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
2131                                                  , AliTRDseedV1/*AliRieman*/ *reference)
2132 {
2133   //
2134   // Finds a seeding Cluster for the extrapolation chamber.
2135   //
2136   // The seeding cluster should be as close as possible to the assumed
2137   // track which is represented by a Rieman fit.
2138   // Therefore the selecting criterion is the minimum distance between
2139   // the best fitting cluster and the Reference which is derived from
2140   // the AliTRDseed. Because all layers are assumed to be equally good
2141   // a linear search is performed.
2142   //
2143   // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!)
2144   //                   - sfit: the reference
2145   //
2146   // Output:           - the best seeding cluster
2147   //
2148
2149         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2150         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2151
2152         // distances as squared distances
2153         Int_t index = 0;
2154         Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearestDistance =100000.0;
2155         ypos = reference->GetYref(0);
2156         zpos = reference->GetZref(0);
2157         AliTRDcluster *currentBest = 0x0, *temp = 0x0;
2158         for(Int_t ils = 0; ils < nTimeBins; ils++){
2159                 // Reference positions
2160 //              ypos = reference->GetYat(layers[ils].GetX());
2161 //              zpos = reference->GetZat(layers[ils].GetX());
2162                 index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z());
2163                 if(index == -1) continue;
2164                 temp = layers[ils].GetCluster(index);
2165                 if(!temp) continue;
2166                 distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
2167                 if(distance < nearestDistance){
2168                         nearestDistance = distance;
2169                         currentBest = temp;
2170                 }
2171         }
2172         return currentBest;
2173 }
2174
2175 //____________________________________________________________________
2176 AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
2177                                                   , Int_t plane)
2178 {
2179   //
2180   // Creates a seeding layer
2181   //
2182
2183         // constants
2184         const Int_t kMaxRows = 16;
2185         const Int_t kMaxCols = 144;
2186         const Int_t kMaxPads = 2304;
2187
2188         // Get the calculation
2189         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
2190         Int_t nTimeBins = cal->GetNumberOfTimeBins();
2191
2192         // Get the geometrical data of the chamber
2193         AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr());
2194         Int_t nCols = pp->GetNcols();
2195         Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
2196         Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
2197         Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
2198         Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
2199         Int_t nRows = pp->GetNrows();
2200         Float_t binlength = (ymax - ymin)/nCols;
2201         //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
2202
2203         // Fill the histogram
2204         Int_t arrpos;
2205         Float_t ypos;
2206         Int_t irow, nClusters;
2207         Int_t *histogram[kMaxRows];                                                                                     // 2D-Histogram
2208         Int_t hvals[kMaxPads];  memset(hvals, 0, sizeof(Int_t)*kMaxPads);
2209         Float_t *sigmas[kMaxRows];
2210         Float_t svals[kMaxPads];        memset(svals, 0, sizeof(Float_t)*kMaxPads);
2211         AliTRDcluster *c = 0x0;
2212         for(Int_t irs = 0; irs < kMaxRows; irs++){
2213                 histogram[irs] = &hvals[irs*kMaxCols];
2214                 sigmas[irs] = &svals[irs*kMaxCols];
2215         }
2216         for(Int_t iTime = 0; iTime < nTimeBins; iTime++){
2217                 nClusters = layers[iTime].GetNClusters();
2218                 for(Int_t incl = 0; incl < nClusters; incl++){
2219                         c = layers[iTime].GetCluster(incl);
2220                         ypos = c->GetY();
2221                         if(ypos > ymax && ypos < ymin) continue;
2222                         irow = pp->GetPadRowNumber(c->GetZ());                          // Zbin
2223                         if(irow < 0)continue;
2224                         arrpos = static_cast<Int_t>((ypos - ymin)/binlength);
2225                         if(ypos == ymax) arrpos = nCols - 1;
2226                         histogram[irow][arrpos]++;
2227                         sigmas[irow][arrpos] += c->GetSigmaZ2();
2228                 }
2229         }
2230
2231 // Now I have everything in the histogram, do the selection
2232 //      printf("Starting the analysis\n");
2233         //Int_t nPads = nCols * nRows;
2234         // This is what we are interested in: The center of gravity of the best candidates
2235         Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
2236         Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
2237         Float_t *cogy[kMaxRows];
2238         Float_t *cogz[kMaxRows];
2239         // Lookup-Table storing coordinates according ti the bins
2240         Float_t yLengths[kMaxCols];
2241         Float_t zLengths[kMaxRows];
2242         for(Int_t icnt = 0; icnt < nCols; icnt++){
2243                 yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
2244         }
2245         for(Int_t icnt = 0; icnt < nRows; icnt++){
2246                 zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
2247         }
2248
2249         // A bitfield is used to mask the pads as usable
2250         Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
2251         for(UChar_t icount = 0; icount < nRows; icount++){
2252                 cogy[icount] = &cogyvals[icount*kMaxCols];
2253                 cogz[icount] = &cogzvals[icount*kMaxCols];
2254         }
2255         // In this array the array position of the best candidates will be stored
2256         Int_t cand[kMaxTracksStack];
2257         Float_t sigcands[kMaxTracksStack];
2258
2259         // helper variables
2260         Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads);
2261         Int_t nCandidates = 0;
2262         Float_t norm, cogv;
2263         // histogram filled -> Select best bins
2264         TMath::Sort(kMaxPads, hvals, indices);                  // bins storing a 0 should not matter
2265         // Set Threshold
2266         Int_t maximum = hvals[indices[0]];      // best
2267         Int_t threshold = static_cast<UChar_t>(maximum * fRecoParam->GetFindableClusters());
2268         Int_t col, row, lower, lower1, upper, upper1;
2269         for(Int_t ib = 0; ib < kMaxPads; ib++){
2270                 if(nCandidates >= kMaxTracksStack){
2271                         AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, kMaxTracksStack));
2272                         break;
2273                 }
2274                 // Positions
2275                 row = indices[ib]/nCols;
2276                 col = indices[ib]%nCols;
2277                 // here will be the threshold condition:
2278                 if((mask[col] & (1 << row)) != 0) continue;             // Pad is masked: continue
2279                 if(histogram[row][col] < TMath::Max(threshold, 1)){     // of course at least one cluster is needed
2280                         break;                  // number of clusters below threshold: break;
2281                 }
2282                 // passing: Mark the neighbors
2283                 lower  = TMath::Max(col - 1, 0); upper  = TMath::Min(col + 2, nCols);
2284                 lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
2285                 for(Int_t ic = lower; ic < upper; ++ic)
2286                         for(Int_t ir = lower1; ir < upper1; ++ir){
2287                                 if(ic == col && ir == row) continue;
2288                                 mask[ic] |= (1 << ir);
2289                         }
2290                 // Storing the position in an array
2291                 // testing for neigboring
2292                 cogv = 0;
2293                 norm = 0;
2294                 lower = TMath::Max(col - 1,0);
2295                 upper = TMath::Min(col + 2, nCols);
2296                 for(Int_t inb = lower; inb < upper; ++inb){
2297                         cogv += yLengths[inb] * histogram[row][inb];
2298                         norm += histogram[row][inb];
2299                 }
2300                 cogy[row][col] = cogv / norm;
2301                 cogv = 0; norm = 0;
2302                 lower = TMath::Max(row - 1, 0);
2303                 upper = TMath::Min(row + 2, nRows);
2304                 for(Int_t inb = lower; inb < upper; ++inb){
2305                         cogv += zLengths[inb] * histogram[inb][col];
2306                         norm += histogram[inb][col];
2307                 }
2308                 cogz[row][col] = cogv /  norm;
2309                 // passed the filter
2310                 cand[nCandidates] = row*kMaxCols + col; // store the position of a passig candidate into an Array
2311                 sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
2312                 // Analysis output
2313                 nCandidates++;
2314         }
2315         AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr());
2316         fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2);
2317         fakeLayer->SetSector(layers[0].GetSector());
2318         AliTRDcluster **fakeClusters = 0x0;
2319         UInt_t *fakeIndices = 0x0;
2320         if(nCandidates){
2321                 fakeClusters = new AliTRDcluster*[nCandidates];
2322                 fakeIndices = new UInt_t[nCandidates];
2323                 UInt_t fakeIndex = 0;
2324                 for(Int_t ican = 0; ican < nCandidates; ican++){
2325                         fakeClusters[ican] = new AliTRDcluster();
2326                         fakeClusters[ican]->SetX(fakeLayer->GetX());
2327                         fakeClusters[ican]->SetY(cogyvals[cand[ican]]);
2328                         fakeClusters[ican]->SetZ(cogzvals[cand[ican]]);
2329                         fakeClusters[ican]->SetSigmaZ2(sigcands[ican]);
2330                         fakeIndices[ican] = fakeIndex++;// fantasy number
2331                 }
2332         }
2333         fakeLayer->SetRecoParam(fRecoParam);
2334         fakeLayer->SetClustersArray(fakeClusters, nCandidates);
2335         fakeLayer->SetIndexArray(fakeIndices);
2336         fakeLayer->SetNRows(nRows);
2337         fakeLayer->BuildIndices();
2338         //fakeLayer->PrintClusters();
2339
2340 #ifdef DEBUG
2341         if(AliTRDReconstructor::StreamLevel() >= 3){
2342                 TMatrixD hist(nRows, nCols);
2343                 for(Int_t i = 0; i < nRows; i++)
2344                         for(Int_t j = 0; j < nCols; j++)
2345                                 hist(i,j) = histogram[i][j];
2346                 TTreeSRedirector &cstreamer = *fDebugStreamer;
2347                 cstreamer << "MakeSeedingLayer"
2348                         << "Iteration="  << fSieveSeeding
2349                         << "plane="      << plane
2350                         << "ymin="       << ymin
2351                         << "ymax="       << ymax
2352                         << "zmin="       << zmin
2353                         << "zmax="       << zmax
2354                         << "L.="         << fakeLayer
2355                         << "Histogram.=" << &hist
2356                         << "\n";
2357         }
2358 #endif
2359         return fakeLayer;
2360 }
2361
2362 //____________________________________________________________________
2363 void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4])
2364 {
2365   //
2366   // Map seeding configurations to detector planes.
2367   //
2368   // Parameters :
2369   //   iconfig : configuration index
2370   //   planes  : member planes of this configuration. On input empty.
2371   //
2372   // Output :
2373   //   planes : contains the planes which are defining the configuration
2374   //
2375   // Detailed description
2376   //
2377   // Here is the list of seeding planes configurations together with
2378   // their topological classification:
2379   //
2380   //  0 - 5432 TQ 0
2381   //  1 - 4321 TQ 0
2382   //  2 - 3210 TQ 0
2383   //  3 - 5321 TQ 1
2384   //  4 - 4210 TQ 1
2385   //  5 - 5431 TQ 1
2386   //  6 - 4320 TQ 1
2387   //  7 - 5430 TQ 2
2388   //  8 - 5210 TQ 2
2389   //  9 - 5421 TQ 3
2390   // 10 - 4310 TQ 3
2391   // 11 - 5410 TQ 4
2392   // 12 - 5420 TQ 5
2393   // 13 - 5320 TQ 5
2394   // 14 - 5310 TQ 5
2395   //
2396   // The topologic quality is modeled as follows:
2397   // 1. The general model is define by the equation:
2398   //  p(conf) = exp(-conf/2)
2399   // 2. According to the topologic classification, configurations from the same
2400   //    class are assigned the agerage value over the model values.
2401   // 3. Quality values are normalized.
2402   //
2403   // The topologic quality distribution as function of configuration is given below:
2404   //Begin_Html
2405   // <img src="gif/topologicQA.gif">
2406   //End_Html
2407   //
2408
2409         switch(iconfig){
2410         case 0: // 5432 TQ 0
2411                 planes[0] = 2;
2412                 planes[1] = 3;
2413                 planes[2] = 4;
2414                 planes[3] = 5;
2415                 break;
2416         case 1: // 4321 TQ 0
2417                 planes[0] = 1;
2418                 planes[1] = 2;
2419                 planes[2] = 3;
2420                 planes[3] = 4;
2421                 break;
2422         case 2: // 3210 TQ 0
2423                 planes[0] = 0;
2424                 planes[1] = 1;
2425                 planes[2] = 2;
2426                 planes[3] = 3;
2427                 break;
2428         case 3: // 5321 TQ 1
2429                 planes[0] = 1;
2430                 planes[1] = 2;
2431                 planes[2] = 3;
2432                 planes[3] = 5;
2433                 break;
2434         case 4: // 4210 TQ 1
2435                 planes[0] = 0;
2436                 planes[1] = 1;
2437                 planes[2] = 2;
2438                 planes[3] = 4;
2439                 break;
2440         case 5: // 5431 TQ 1
2441                 planes[0] = 1;
2442                 planes[1] = 3;
2443                 planes[2] = 4;
2444                 planes[3] = 5;
2445                 break;
2446         case 6: // 4320 TQ 1
2447                 planes[0] = 0;
2448                 planes[1] = 2;
2449                 planes[2] = 3;
2450                 planes[3] = 4;
2451                 break;
2452         case 7: // 5430 TQ 2
2453                 planes[0] = 0;
2454                 planes[1] = 3;
2455                 planes[2] = 4;
2456                 planes[3] = 5;
2457                 break;
2458         case 8: // 5210 TQ 2
2459                 planes[0] = 0;
2460                 planes[1] = 1;
2461                 planes[2] = 2;
2462                 planes[3] = 5;
2463                 break;
2464         case 9: // 5421 TQ 3
2465                 planes[0] = 1;
2466                 planes[1] = 2;
2467                 planes[2] = 4;
2468                 planes[3] = 5;
2469                 break;
2470         case 10: // 4310 TQ 3
2471                 planes[0] = 0;
2472                 planes[1] = 1;
2473                 planes[2] = 3;
2474                 planes[3] = 4;
2475                 break;
2476         case 11: // 5410 TQ 4
2477                 planes[0] = 0;
2478                 planes[1] = 1;
2479                 planes[2] = 4;
2480                 planes[3] = 5;
2481                 break;
2482         case 12: // 5420 TQ 5
2483                 planes[0] = 0;
2484                 planes[1] = 2;
2485                 planes[2] = 4;
2486                 planes[3] = 5;
2487                 break;
2488         case 13: // 5320 TQ 5
2489                 planes[0] = 0;
2490                 planes[1] = 2;
2491                 planes[2] = 3;
2492                 planes[3] = 5;
2493                 break;
2494         case 14: // 5310 TQ 5
2495                 planes[0] = 0;
2496                 planes[1] = 1;
2497                 planes[2] = 3;
2498                 planes[3] = 5;
2499                 break;
2500         }
2501 }
2502
2503 //____________________________________________________________________
2504 void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
2505 {
2506   //
2507   // Returns the extrapolation planes for a seeding configuration.
2508   //
2509   // Parameters :
2510   //   iconfig : configuration index
2511   //   planes  : planes which are not in this configuration. On input empty.
2512   //
2513   // Output :
2514   //   planes : contains the planes which are not in the configuration
2515   //
2516   // Detailed description
2517   //
2518
2519         switch(iconfig){
2520         case 0: // 5432 TQ 0
2521                 planes[0] = 1;
2522                 planes[1] = 0;
2523                 break;
2524         case 1: // 4321 TQ 0
2525                 planes[0] = 5;
2526                 planes[1] = 0;
2527                 break;
2528         case 2: // 3210 TQ 0
2529                 planes[0] = 4;
2530                 planes[1] = 5;
2531                 break;
2532         case 3: // 5321 TQ 1
2533                 planes[0] = 4;
2534                 planes[1] = 0;
2535                 break;
2536         case 4: // 4210 TQ 1
2537                 planes[0] = 5;
2538                 planes[1] = 3;
2539                 break;
2540         case 5: // 5431 TQ 1
2541                 planes[0] = 2;
2542                 planes[1] = 0;
2543                 break;
2544         case 6: // 4320 TQ 1
2545                 planes[0] = 5;
2546                 planes[1] = 1;
2547                 break;
2548         case 7: // 5430 TQ 2
2549                 planes[0] = 2;
2550                 planes[1] = 1;
2551                 break;
2552         case 8: // 5210 TQ 2
2553                 planes[0] = 4;
2554                 planes[1] = 3;
2555                 break;
2556         case 9: // 5421 TQ 3
2557                 planes[0] = 3;
2558                 planes[1] = 0;
2559                 break;
2560         case 10: // 4310 TQ 3
2561                 planes[0] = 5;
2562                 planes[1] = 2;
2563                 break;
2564         case 11: // 5410 TQ 4
2565                 planes[0] = 3;
2566                 planes[1] = 2;
2567                 break;
2568         case 12: // 5420 TQ 5
2569                 planes[0] = 3;
2570                 planes[1] = 1;
2571                 break;
2572         case 13: // 5320 TQ 5
2573                 planes[0] = 4;
2574                 planes[1] = 1;
2575                 break;
2576         case 14: // 5310 TQ 5
2577                 planes[0] = 4;
2578                 planes[1] = 2;
2579                 break;
2580         }
2581 }