TRD/AliTRDtrackerV1.cxx

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