TRD/AliTRDclusterizer.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id$ */
  18
  19 ///////////////////////////////////////////////////////////////////////////////
  20 //                                                                           //
  21 // TRD cluster finder                                                        //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include <TF1.h>
  26 #include <TTree.h>
  27 #include <TH1.h>
  28 #include <TFile.h>
  29 #include <TObjArray.h>
  30
  31 #include "AliRunLoader.h"
  32 #include "AliLoader.h"
  33 #include "AliRawReader.h"
  34 #include "AliLog.h"
  35 #include "AliAlignObj.h"
  36
  37 #include "AliTRDclusterizer.h"
  38 #include "AliTRDcluster.h"
  39 #include "AliTRDReconstructor.h"
  40 #include "AliTRDgeometry.h"
  41 #include "AliTRDdataArrayF.h"
  42 #include "AliTRDdataArrayI.h"
  43 #include "AliTRDdataArrayS.h"
  44 #include "AliTRDdigitsManager.h"
  45 #include "AliTRDrawData.h"
  46 #include "AliTRDcalibDB.h"
  47 #include "AliTRDrecoParam.h"
  48 #include "AliTRDCommonParam.h"
  49 #include "AliTRDtransform.h"
  50 #include "AliTRDSignalIndex.h"
  51 #include "AliTRDrawStreamBase.h"
  52 #include "AliTRDfeeParam.h"
  53
  54 #include "Cal/AliTRDCalROC.h"
  55 #include "Cal/AliTRDCalDet.h"
  56
  57 ClassImp(AliTRDclusterizer)
  58
  59 //_____________________________________________________________________________
  60 AliTRDclusterizer::AliTRDclusterizer()
  61   :TNamed()
  62   ,fRunLoader(NULL)
  63   ,fClusterTree(NULL)
  64   ,fRecPoints(NULL)
  65   ,fDigitsManager(NULL)
  66   ,fAddLabels(kTRUE)
  67   ,fRawVersion(2)
  68   ,fIndexesOut(NULL)
  69   ,fIndexesMaxima(NULL)
  70   ,fTransform(new AliTRDtransform(0))
  71   ,fLUTbin(0)
  72   ,fLUT(NULL)
  73 {
  74   //
  75   // AliTRDclusterizer default constructor
  76   //
  77
  78   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
  79
  80 }
  81
  82 //_____________________________________________________________________________
  83 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title)
  84   :TNamed(name,title)
  85   ,fRunLoader(NULL)
  86   ,fClusterTree(NULL)
  87   ,fRecPoints(NULL)
  88   ,fDigitsManager(new AliTRDdigitsManager())
  89   ,fAddLabels(kTRUE)
  90   ,fRawVersion(2)
  91   ,fIndexesOut(NULL)
  92   ,fIndexesMaxima(NULL)
  93   ,fTransform(new AliTRDtransform(0))
  94   ,fLUTbin(0)
  95   ,fLUT(NULL)
  96 {
  97   //
  98   // AliTRDclusterizer constructor
  99   //
 100
 101   fDigitsManager->CreateArrays();
 102
 103   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 104
 105   FillLUT();
 106
 107 }
 108
 109 //_____________________________________________________________________________
 110 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
 111   :TNamed(c)
 112   ,fRunLoader(NULL)
 113   ,fClusterTree(NULL)
 114   ,fRecPoints(NULL)
 115   ,fDigitsManager(NULL)
 116   ,fAddLabels(kTRUE)
 117   ,fRawVersion(2)
 118   ,fIndexesOut(NULL)
 119   ,fIndexesMaxima(NULL)
 120   ,fTransform(NULL)
 121   ,fLUTbin(0)
 122   ,fLUT(0)
 123 {
 124   //
 125   // AliTRDclusterizer copy constructor
 126   //
 127
 128   FillLUT();
 129
 130 }
 131
 132 //_____________________________________________________________________________
 133 AliTRDclusterizer::~AliTRDclusterizer()
 134 {
 135   //
 136   // AliTRDclusterizer destructor
 137   //
 138
 139   if (fRecPoints)
 140     {
 141       fRecPoints->Delete();
 142       delete fRecPoints;
 143     }
 144
 145   if (fDigitsManager)
 146     {
 147       delete fDigitsManager;
 148       fDigitsManager = NULL;
 149     }
 150
 151   if (fIndexesOut)
 152     {
 153       delete fIndexesOut;
 154       fIndexesOut    = NULL;
 155     }
 156
 157   if (fIndexesMaxima)
 158     {
 159       delete fIndexesMaxima;
 160       fIndexesMaxima = NULL;
 161     }
 162
 163   if (fTransform)
 164     {
 165       delete fTransform;
 166       fTransform     = NULL;
 167     }
 168
 169   if (fLUT)
 170     {
 171       delete [] fLUT;
 172       fLUT           = NULL;
 173     }
 174
 175 }
 176
 177 //_____________________________________________________________________________
 178 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
 179 {
 180   //
 181   // Assignment operator
 182   //
 183
 184   if (this != &c)
 185     {
 186       ((AliTRDclusterizer &) c).Copy(*this);
 187     }
 188   return *this;
 189
 190 }
 191
 192 //_____________________________________________________________________________
 193 void AliTRDclusterizer::Copy(TObject &c) const
 194 {
 195   //
 196   // Copy function
 197   //
 198
 199   ((AliTRDclusterizer &) c).fClusterTree   = NULL;
 200   ((AliTRDclusterizer &) c).fRecPoints     = NULL;
 201   ((AliTRDclusterizer &) c).fDigitsManager = NULL;
 202   ((AliTRDclusterizer &) c).fAddLabels     = fAddLabels;
 203   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
 204   ((AliTRDclusterizer &) c).fIndexesOut    = NULL;
 205   ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
 206   ((AliTRDclusterizer &) c).fTransform     = NULL;
 207   ((AliTRDclusterizer &) c).fLUTbin        = 0;
 208   ((AliTRDclusterizer &) c).fLUT           = NULL;
 209
 210 }
 211
 212 //_____________________________________________________________________________
 213 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
 214 {
 215   //
 216   // Opens the AliROOT file. Output and input are in the same file
 217   //
 218
 219   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 220   fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 221
 222   if (!fRunLoader) {
 223     fRunLoader = AliRunLoader::Open(name);
 224   }
 225
 226   if (!fRunLoader) {
 227     AliError(Form("Can not open session for file %s.",name));
 228     return kFALSE;
 229   }
 230
 231   OpenInput(nEvent);
 232   OpenOutput();
 233
 234   return kTRUE;
 235
 236 }
 237
 238 //_____________________________________________________________________________
 239 Bool_t AliTRDclusterizer::OpenOutput()
 240 {
 241   //
 242   // Open the output file
 243   //
 244
 245   TObjArray *ioArray = 0;
 246
 247   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 248   loader->MakeTree("R");
 249
 250   fClusterTree = loader->TreeR();
 251   fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
 252
 253   return kTRUE;
 254
 255 }
 256
 257 //_____________________________________________________________________________
 258 Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
 259 {
 260   //
 261   // Connect the output tree
 262   //
 263
 264   TObjArray *ioArray = 0;
 265
 266   fClusterTree = clusterTree;
 267   fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
 268
 269   return kTRUE;
 270
 271 }
 272
 273 //_____________________________________________________________________________
 274 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
 275 {
 276   //
 277   // Opens a ROOT-file with TRD-hits and reads in the digits-tree
 278   //
 279
 280   // Import the Trees for the event nEvent in the file
 281   fRunLoader->GetEvent(nEvent);
 282
 283   return kTRUE;
 284
 285 }
 286
 287 //_____________________________________________________________________________
 288 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
 289 {
 290   //
 291   // Fills TRDcluster branch in the tree with the clusters
 292   // found in detector = det. For det=-1 writes the tree.
 293   //
 294
 295   if ((det <                      -1) ||
 296       (det >= AliTRDgeometry::Ndet())) {
 297     AliError(Form("Unexpected detector index %d.\n",det));
 298     return kFALSE;
 299   }
 300
 301   TBranch *branch = fClusterTree->GetBranch("TRDcluster");
 302   if (!branch) {
 303     TObjArray *ioArray = 0;
 304     branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
 305   }
 306
 307   if ((det >=                      0) &&
 308       (det <  AliTRDgeometry::Ndet())) {
 309
 310     Int_t nRecPoints = RecPoints()->GetEntriesFast();
 311     TObjArray *detRecPoints = new TObjArray(400);
 312
 313     for (Int_t i = 0; i < nRecPoints; i++) {
 314       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 315       if (det == c->GetDetector()) {
 316         detRecPoints->AddLast(c);
 317       }
 318       else {
 319         AliError(Form("Attempt to write a cluster with unexpected detector index: got=%d expected=%d\n"
 320                      ,c->GetDetector()
 321                      ,det));
 322       }
 323     }
 324
 325     branch->SetAddress(&detRecPoints);
 326     fClusterTree->Fill();
 327
 328     delete detRecPoints;
 329
 330     return kTRUE;
 331
 332   }
 333
 334   if (det == -1) {
 335
 336     AliInfo(Form("Writing the cluster tree %s for event %d."
 337                 ,fClusterTree->GetName(),fRunLoader->GetEventNumber()));
 338
 339     if (fRecPoints) {
 340
 341       branch->SetAddress(&fRecPoints);
 342
 343       AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
 344       loader->WriteRecPoints("OVERWRITE");
 345
 346     }
 347     else {
 348
 349       AliError("Cluster tree does not exist. Cannot write clusters.\n");
 350       return kFALSE;
 351
 352     }
 353
 354     return kTRUE;
 355
 356   }
 357
 358   AliError(Form("Unexpected detector index %d.\n",det));
 359
 360   return kFALSE;
 361
 362 }
 363
 364 //_____________________________________________________________________________
 365 void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
 366 {
 367   //
 368   // Reset the helper indexes
 369   //
 370
 371   if (fIndexesOut)
 372     {
 373       // carefull here - we assume that only row number may change - most probable
 374       if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
 375         fIndexesOut->ResetContent();
 376       else
 377         fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
 378                                            , indexesIn->GetNcol()
 379                                            , indexesIn->GetNtime());
 380     }
 381   else
 382     {
 383       fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
 384                                         , indexesIn->GetNcol()
 385                                         , indexesIn->GetNtime());
 386     }
 387
 388   if (fIndexesMaxima)
 389     {
 390       // carefull here - we assume that only row number may change - most probable
 391       if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
 392         {
 393           fIndexesMaxima->ResetContent();
 394         }
 395       else
 396         {
 397           fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
 398                                                 , indexesIn->GetNcol()
 399                                                 , indexesIn->GetNtime());
 400         }
 401     }
 402   else
 403     {
 404       fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
 405                                            , indexesIn->GetNcol()
 406                                            , indexesIn->GetNtime());
 407     }
 408
 409 }
 410
 411 //_____________________________________________________________________________
 412 Bool_t AliTRDclusterizer::ReadDigits()
 413 {
 414   //
 415   // Reads the digits arrays from the input aliroot file
 416   //
 417
 418   if (!fRunLoader) {
 419     AliError("No run loader available");
 420     return kFALSE;
 421   }
 422
 423   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 424   if (!loader->TreeD()) {
 425     loader->LoadDigits();
 426   }
 427
 428   // Read in the digit arrays
 429   return (fDigitsManager->ReadDigits(loader->TreeD()));
 430
 431 }
 432
 433 //_____________________________________________________________________________
 434 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
 435 {
 436   //
 437   // Reads the digits arrays from the input tree
 438   //
 439
 440   // Read in the digit arrays
 441   return (fDigitsManager->ReadDigits(digitsTree));
 442
 443 }
 444
 445 //_____________________________________________________________________________
 446 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
 447 {
 448   //
 449   // Reads the digits arrays from the ddl file
 450   //
 451
 452   AliTRDrawData raw;
 453   fDigitsManager = raw.Raw2Digits(rawReader);
 454
 455   return kTRUE;
 456
 457 }
 458
 459 //_____________________________________________________________________________
 460 Bool_t AliTRDclusterizer::MakeClusters()
 461 {
 462   //
 463   // Creates clusters from digits
 464   //
 465
 466   // Propagate info from the digits manager
 467   if (fAddLabels == kTRUE)
 468     {
 469       fAddLabels = fDigitsManager->UsesDictionaries();
 470     }
 471
 472   Bool_t fReturn = kTRUE;
 473   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
 474     {
 475
 476       AliTRDdataArrayS *digitsIn = (AliTRDdataArrayS *) fDigitsManager->GetDigits(i);
 477       // This is to take care of switched off super modules
 478       if (!digitsIn->HasData())
 479         {
 480           continue;
 481         }
 482       digitsIn->Expand();
 483       AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 484       if (indexes->IsAllocated() == kFALSE)
 485         {
 486           fDigitsManager->BuildIndexes(i);
 487         }
 488
 489       Bool_t fR = kFALSE;
 490       if (indexes->HasEntry())
 491         {
 492           if (fAddLabels)
 493             {
 494               for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
 495                 {
 496                   AliTRDdataArrayI *tracksIn = 0;
 497                   tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(i,iDict);
 498                   tracksIn->Expand();
 499                 }
 500             }
 501           fR = MakeClusters(i);
 502           fReturn = fR && fReturn;
 503         }
 504
 505       if (fR == kFALSE)
 506         {
 507           WriteClusters(i);
 508           ResetRecPoints();
 509         }
 510
 511       // No compress just remove
 512       fDigitsManager->RemoveDigits(i);
 513       fDigitsManager->RemoveDictionaries(i);
 514       fDigitsManager->ClearIndexes(i);
 515
 516     }
 517
 518   return fReturn;
 519
 520 }
 521
 522 //_____________________________________________________________________________
 523 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
 524 {
 525   //
 526   // Creates clusters from raw data
 527   //
 528
 529   return Raw2ClustersChamber(rawReader);
 530
 531 }
 532
 533 //_____________________________________________________________________________
 534 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 535 {
 536   //
 537   // Creates clusters from raw data
 538   //
 539
 540   // Create the digits manager
 541   if (!fDigitsManager)
 542     {
 543       fDigitsManager = new AliTRDdigitsManager();
 544       fDigitsManager->CreateArrays();
 545     }
 546
 547   fDigitsManager->SetUseDictionaries(fAddLabels);
 548
 549   AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
 550   AliTRDrawStreamBase &input = *pinput;
 551
 552   AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
 553
 554   Int_t det    = 0;
 555   while ((det = input.NextChamber(fDigitsManager)) >= 0)
 556     {
 557       Bool_t iclusterBranch = kFALSE;
 558       if (fDigitsManager->GetIndexes(det)->HasEntry())
 559         {
 560           iclusterBranch = MakeClusters(det);
 561         }
 562       if (iclusterBranch == kFALSE)
 563         {
 564           WriteClusters(det);
 565           ResetRecPoints();
 566         }
 567       fDigitsManager->RemoveDigits(det);
 568       fDigitsManager->RemoveDictionaries(det);
 569       fDigitsManager->ClearIndexes(det);
 570     }
 571
 572   delete fDigitsManager;
 573   fDigitsManager = NULL;
 574
 575   delete pinput;
 576   pinput = NULL;
 577   return kTRUE;
 578
 579 }
 580
 581 //_____________________________________________________________________________
 582 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 583 {
 584   //
 585   // Generates the cluster.
 586   //
 587
 588   // Get the digits
 589   //   digits should be expanded beforehand!
 590   //   digitsIn->Expand();
 591   AliTRDdataArrayS *digitsIn = (AliTRDdataArrayS *) fDigitsManager->GetDigits(det);
 592
 593   // This is to take care of switched off super modules
 594   if (!digitsIn->HasData())
 595     {
 596       return kFALSE;
 597     }
 598
 599   AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
 600   if (indexesIn->IsAllocated() == kFALSE)
 601     {
 602       AliError("Indexes do not exist!");
 603       return kFALSE;
 604     }
 605
 606   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
 607   if (!calibration)
 608     {
 609       AliFatal("No AliTRDcalibDB instance available\n");
 610       return kFALSE;
 611     }
 612
 613   // ADC thresholds
 614   // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
 615   Float_t adcThreshold   = 0;
 616
 617   if (!AliTRDReconstructor::RecoParam())
 618     {
 619       AliError("RecoParam does not exist\n");
 620       return kFALSE;
 621     }
 622
 623   // Threshold value for the maximum
 624   Float_t maxThresh      = AliTRDReconstructor::RecoParam()->GetClusMaxThresh();
 625   // Threshold value for the digit signal
 626   Float_t sigThresh      = AliTRDReconstructor::RecoParam()->GetClusSigThresh();
 627
 628   // Iteration limit for unfolding procedure
 629   const Float_t kEpsilon = 0.01;
 630   const Int_t   kNclus   = 3;
 631   const Int_t   kNsig    = 5;
 632
 633   Int_t    iUnfold       = 0;
 634   Double_t ratioLeft     = 1.0;
 635   Double_t ratioRight    = 1.0;
 636
 637   Double_t padSignal[kNsig];
 638   Double_t clusterSignal[kNclus];
 639
 640   Int_t icham = indexesIn->GetChamber();
 641   Int_t iplan = indexesIn->GetPlane();
 642   Int_t isect = indexesIn->GetSM();
 643
 644   // Start clustering in the chamber
 645
 646   Int_t idet  = AliTRDgeometry::GetDetector(iplan,icham,isect);
 647   if (idet != det)
 648     {
 649       AliError("Strange Detector number Missmatch!");
 650       return kFALSE;
 651     }
 652
 653   // TRD space point transformation
 654   fTransform->SetDetector(det);
 655
 656   Int_t    ilayer  = AliGeomManager::kTRD1 + iplan;
 657   Int_t    imodule = icham + AliTRDgeometry::Ncham() * isect;
 658   UShort_t volid   = AliGeomManager::LayerToVolUID(ilayer,imodule);
 659
 660   Int_t nColMax    = digitsIn->GetNcol();
 661   Int_t nTimeTotal = digitsIn->GetNtime();
 662
 663   // Detector wise calibration object for the gain factors
 664   const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();
 665   // Calibration object with pad wise values for the gain factors
 666   AliTRDCalROC       *calGainFactorROC      = calibration->GetGainFactorROC(idet);
 667   // Calibration value for chamber wise gain factor
 668   Float_t             calGainFactorDetValue = calGainFactorDet->GetValue(idet);
 669
 670   Int_t nClusters = 0;
 671
 672   AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow()
 673                                                     ,digitsIn->GetNcol()
 674                                                     ,digitsIn->GetNtime());
 675
 676   ResetHelperIndexes(indexesIn);
 677
 678   // Apply the gain and the tail cancelation via digital filter
 679   TailCancelation(digitsIn
 680                  ,digitsOut
 681                  ,indexesIn
 682                  ,fIndexesOut
 683                  ,nTimeTotal
 684                  ,adcThreshold
 685                  ,calGainFactorROC
 686                  ,calGainFactorDetValue);
 687
 688   Int_t row  = 0;
 689   Int_t col  = 0;
 690   Int_t time = 0;
 691   Int_t iPad = 0;
 692
 693   fIndexesOut->ResetCounters();
 694   while (fIndexesOut->NextRCTbinIndex(row, col, time))
 695     {
 696
 697       Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
 698
 699       // Look for the maximum
 700       if (signalM >= maxThresh)
 701         {
 702
 703           if (col + 1 >= nColMax || col-1 < 0)
 704             continue;
 705
 706           Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
 707           Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
 708
 709           if ((TMath::Abs(signalL) <= signalM) &&
 710               (TMath::Abs(signalR) <  signalM))
 711             {
 712               if ((TMath::Abs(signalL) >= sigThresh) ||
 713                   (TMath::Abs(signalR) >= sigThresh))
 714                 {
 715                   // Maximum found, mark the position by a negative signal
 716                   digitsOut->SetDataUnchecked(row,col,time,-signalM);
 717                   fIndexesMaxima->AddIndexTBin(row,col,time);
 718                 }
 719             }
 720
 721         }
 722
 723     }
 724
 725   // The index to the first cluster of a given ROC
 726   Int_t firstClusterROC = -1;
 727   // The number of cluster in a given ROC
 728   Int_t nClusterROC     =  0;
 729
 730   // Now check the maxima and calculate the cluster position
 731   fIndexesMaxima->ResetCounters();
 732   while (fIndexesMaxima->NextRCTbinIndex(row, col, time))
 733     {
 734
 735       // Maximum found ?
 736       if (digitsOut->GetDataUnchecked(row,col,time) < 0.0)
 737         {
 738
 739           for (iPad = 0; iPad < kNclus; iPad++)
 740             {
 741               Int_t iPadCol = col - 1 + iPad;
 742               clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
 743             }
 744
 745           // Count the number of pads in the cluster
 746           Int_t nPadCount = 0;
 747           Int_t ii;
 748           // Look to the left
 749           ii = 0;
 750           while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii  ,time)) >= sigThresh)
 751             {
 752               nPadCount++;
 753               ii++;
 754               if (col-ii   <        0) break;
 755             }
 756           // Look to the right
 757           ii = 0;
 758           while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh)
 759             {
 760               nPadCount++;
 761               ii++;
 762               if (col+ii+1 >= nColMax) break;
 763             }
 764           nClusters++;
 765
 766           // Look for 5 pad cluster with minimum in the middle
 767           Bool_t fivePadCluster = kFALSE;
 768           if (col < (nColMax - 3))
 769             {
 770               if (digitsOut->GetDataUnchecked(row,col+2,time) < 0)
 771                 {
 772                   fivePadCluster = kTRUE;
 773                 }
 774               if ((fivePadCluster) && (col < (nColMax - 5)))
 775                 {
 776                   if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh)
 777                     {
 778                       fivePadCluster = kFALSE;
 779                     }
 780                 }
 781               if ((fivePadCluster) && (col >             1))
 782                 {
 783                   if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh)
 784                     {
 785                       fivePadCluster = kFALSE;
 786                     }
 787                 }
 788             }
 789
 790           // 5 pad cluster
 791           // Modify the signal of the overlapping pad for the left part
 792           // of the cluster which remains from a previous unfolding
 793           if (iUnfold)
 794             {
 795               clusterSignal[0] *= ratioLeft;
 796               iUnfold = 0;
 797             }
 798
 799           // Unfold the 5 pad cluster
 800           if (fivePadCluster)
 801             {
 802               for (iPad = 0; iPad < kNsig; iPad++)
 803                 {
 804                   padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
 805                                                                           ,col-1+iPad
 806                                                                           ,time));
 807                 }
 808               // Unfold the two maxima and set the signal on
 809               // the overlapping pad to the ratio
 810               ratioRight        = Unfold(kEpsilon,iplan,padSignal);
 811               ratioLeft         = 1.0 - ratioRight;
 812               clusterSignal[2] *= ratioRight;
 813               iUnfold = 1;
 814             }
 815
 816           // The position of the cluster in COL direction relative to the center pad (pad units)
 817           Double_t clusterPosCol = 0.0;
 818           if (AliTRDReconstructor::RecoParam()->LUTOn())
 819             {
 820               // Calculate the position of the cluster by using the
 821               // lookup table method
 822               clusterPosCol = LUTposition(iplan,clusterSignal[0]
 823                                                ,clusterSignal[1]
 824                                                ,clusterSignal[2]);
 825             }
 826           else
 827             {
 828               // Calculate the position of the cluster by using the
 829               // center of gravity method
 830               for (Int_t i = 0; i < kNsig; i++)
 831                 {
 832                   padSignal[i] = 0.0;
 833                 }
 834               padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col  ,time)); // Central pad
 835               padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left    pad
 836               padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right   pad
 837               if ((col >           2) &&
 838                   (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1]))
 839                 {
 840                   padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
 841                 }
 842               if ((col < nColMax - 3) &&
 843                   (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3]))
 844                 {
 845                   padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
 846                 }
 847               clusterPosCol = GetCOG(padSignal);
 848             }
 849
 850           // Store the amplitudes of the pads in the cluster for later analysis
 851           Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
 852           for (Int_t jPad = col-3; jPad <= col+3; jPad++)
 853             {
 854               if ((jPad <          0) ||
 855                   (jPad >= nColMax-1))
 856                 {
 857                   continue;
 858                 }
 859               signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
 860             }
 861
 862           // Transform the local cluster coordinates into calibrated
 863           // space point positions defined in the local tracking system.
 864           // Here the calibration for T0, Vdrift and ExB is applied as well.
 865           Double_t clusterXYZ[6];
 866           clusterXYZ[0] = clusterPosCol;
 867           clusterXYZ[1] = clusterSignal[0];
 868           clusterXYZ[2] = clusterSignal[1];
 869           clusterXYZ[3] = clusterSignal[2];
 870           clusterXYZ[4] = 0.0;
 871           clusterXYZ[5] = 0.0;
 872           Int_t    clusterRCT[3];
 873           clusterRCT[0] = row;
 874           clusterRCT[1] = col;
 875           clusterRCT[2] = 0;
 876
 877                 Bool_t out = kTRUE;
 878           if (fTransform->Transform(clusterXYZ, clusterRCT, ((UInt_t) time), out, 0)) {
 879
 880             // Add the cluster to the output array
 881             // The track indices will be stored later
 882             Float_t clusterPos[3];
 883             clusterPos[0] = clusterXYZ[0];
 884             clusterPos[1] = clusterXYZ[1];
 885             clusterPos[2] = clusterXYZ[2];
 886             Float_t clusterSig[2];
 887             clusterSig[0] = clusterXYZ[4];
 888             clusterSig[1] = clusterXYZ[5];
 889             Double_t clusterCharge  = clusterXYZ[3];
 890             Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
 891             AliTRDcluster *cluster = new AliTRDcluster(idet
 892                                                       ,clusterCharge
 893                                                       ,clusterPos
 894                                                       ,clusterSig
 895                                                       ,0x0
 896                                                       ,((Char_t) nPadCount)
 897                                                       ,signals
 898                                                       ,((UChar_t) col)
 899                                                       ,((UChar_t) row)
 900                                                       ,((UChar_t) time)
 901                                                       ,clusterTimeBin
 902                                                       ,clusterPosCol
 903                                                       ,volid);
 904                         cluster->SetInChamber(!out);
 905
 906             // Temporarily store the row, column and time bin of the center pad
 907             // Used to later on assign the track indices
 908             cluster->SetLabel( row,0);
 909             cluster->SetLabel( col,1);
 910             cluster->SetLabel(time,2);
 911
 912             RecPoints()->Add(cluster);
 913
 914             // Store the index of the first cluster in the current ROC
 915             if (firstClusterROC < 0)
 916               {
 917                 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
 918               }
 919
 920             // Count the number of cluster in the current ROC
 921             nClusterROC++;
 922
 923           } // if: Transform ok ?
 924
 925         } // if: Maximum found ?
 926
 927     }
 928
 929   delete digitsOut;
 930
 931   if (fAddLabels)
 932     {
 933       AddLabels(idet, firstClusterROC, nClusterROC);
 934     }
 935
 936   // Write the cluster and reset the array
 937   WriteClusters(idet);
 938   ResetRecPoints();
 939
 940   return kTRUE;
 941
 942 }
 943
 944 //_____________________________________________________________________________
 945 Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
 946 {
 947   //
 948   // Add the track indices to the found clusters
 949   //
 950
 951   const Int_t   kNclus  = 3;
 952   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
 953   const Int_t   kNtrack = kNdict * kNclus;
 954
 955   Int_t iClusterROC = 0;
 956
 957   Int_t row  = 0;
 958   Int_t col  = 0;
 959   Int_t time = 0;
 960   Int_t iPad = 0;
 961
 962   // Temporary array to collect the track indices
 963   Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
 964
 965   // Loop through the dictionary arrays one-by-one
 966   // to keep memory consumption low
 967   AliTRDdataArrayI *tracksIn = 0;
 968   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
 969
 970     // tracksIn should be expanded beforehand!
 971     tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
 972
 973     // Loop though the clusters found in this ROC
 974     for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
 975
 976       AliTRDcluster *cluster = (AliTRDcluster *)
 977         RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
 978       row  = cluster->GetLabel(0);
 979       col  = cluster->GetLabel(1);
 980       time = cluster->GetLabel(2);
 981
 982       for (iPad = 0; iPad < kNclus; iPad++) {
 983         Int_t iPadCol = col - 1 + iPad;
 984         Int_t index   = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
 985         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
 986       }
 987
 988     }
 989
 990   }
 991
 992   // Copy the track indices into the cluster
 993   // Loop though the clusters found in this ROC
 994   for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
 995
 996     AliTRDcluster *cluster = (AliTRDcluster *)
 997       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
 998     cluster->SetLabel(-9999,0);
 999     cluster->SetLabel(-9999,1);
1000     cluster->SetLabel(-9999,2);
1001
1002     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1003
1004   }
1005
1006   delete [] idxTracks;
1007
1008   return kTRUE;
1009
1010 }
1011
1012 //_____________________________________________________________________________
1013 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
1014 {
1015   //
1016   // Get COG position
1017   // Used for clusters with more than 3 pads - where LUT not applicable
1018   //
1019
1020   Double_t sum = signal[0]
1021                + signal[1]
1022                + signal[2]
1023                + signal[3]
1024                + signal[4];
1025
1026   Double_t res = (0.0 * (-signal[0] + signal[4])
1027                       + (-signal[1] + signal[3])) / sum;
1028
1029   return res;
1030
1031 }
1032
1033 //_____________________________________________________________________________
1034 Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
1035 {
1036   //
1037   // Method to unfold neighbouring maxima.
1038   // The charge ratio on the overlapping pad is calculated
1039   // until there is no more change within the range given by eps.
1040   // The resulting ratio is then returned to the calling method.
1041   //
1042
1043   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1044   if (!calibration) {
1045     AliError("No AliTRDcalibDB instance available\n");
1046     return kFALSE;
1047   }
1048
1049   Int_t   irc                = 0;
1050   Int_t   itStep             = 0;                 // Count iteration steps
1051
1052   Double_t ratio             = 0.5;               // Start value for ratio
1053   Double_t prevRatio         = 0.0;               // Store previous ratio
1054
1055   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1056   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1057   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1058
1059   // Start the iteration
1060   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1061
1062     itStep++;
1063     prevRatio = ratio;
1064
1065     // Cluster position according to charge ratio
1066     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1067                       / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1068     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1069                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1070
1071     // Set cluster charge ratio
1072     irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
1073     Double_t ampLeft  = padSignal[1] / newSignal[1];
1074     irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
1075     Double_t ampRight = padSignal[3] / newSignal[1];
1076
1077     // Apply pad response to parameters
1078     irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
1079     irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
1080
1081     // Calculate new overlapping ratio
1082     ratio = TMath::Min((Double_t) 1.0
1083                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1084
1085   }
1086
1087   return ratio;
1088
1089 }
1090
1091 //_____________________________________________________________________________
1092 void AliTRDclusterizer::TailCancelation(AliTRDdataArrayS *digitsIn
1093                                       , AliTRDdataArrayF *digitsOut
1094                                       , AliTRDSignalIndex *indexesIn
1095                                       , AliTRDSignalIndex *indexesOut
1096                                       , Int_t nTimeTotal
1097                                       , Float_t adcThreshold
1098                                       , AliTRDCalROC *calGainFactorROC
1099                                       , Float_t calGainFactorDetValue)
1100 {
1101   //
1102   // Applies the tail cancelation and gain factors:
1103   // Transform digitsIn to digitsOut
1104   //
1105
1106   Int_t iRow  = 0;
1107   Int_t iCol  = 0;
1108   Int_t iTime = 0;
1109
1110   Double_t *inADC  = new Double_t[nTimeTotal];  // ADC data before tail cancellation
1111   Double_t *outADC = new Double_t[nTimeTotal];  // ADC data after tail cancellation
1112   indexesIn->ResetCounters();
1113   while (indexesIn->NextRCIndex(iRow, iCol))
1114     {
1115       Float_t  calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1116       Double_t gain                  = calGainFactorDetValue
1117                                      * calGainFactorROCValue;
1118
1119       for (iTime = 0; iTime < nTimeTotal; iTime++)
1120         {
1121           // Apply gain gain factor
1122           inADC[iTime]   = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1123           inADC[iTime]  /= gain;
1124           outADC[iTime]  = inADC[iTime];
1125         }
1126
1127       // Apply the tail cancelation via the digital filter
1128       if (AliTRDReconstructor::RecoParam()->TCOn())
1129         {
1130           DeConvExp(inADC,outADC,nTimeTotal,AliTRDReconstructor::RecoParam()->GetTCnexp());
1131         }
1132
1133       indexesIn->ResetTbinCounter();
1134       while (indexesIn->NextTbinIndex(iTime))
1135         {
1136           // Store the amplitude of the digit if above threshold
1137           if (outADC[iTime] > adcThreshold)
1138             {
1139               digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1140               indexesOut->AddIndexTBin(iRow,iCol,iTime);
1141             }
1142         } // while itime
1143
1144     } // while irow icol
1145
1146   delete [] inADC;
1147   delete [] outADC;
1148
1149   return;
1150
1151 }
1152
1153 //_____________________________________________________________________________
1154 void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1155                                  , Int_t n, Int_t nexp)
1156 {
1157   //
1158   // Tail cancellation by deconvolution for PASA v4 TRF
1159   //
1160
1161   Double_t rates[2];
1162   Double_t coefficients[2];
1163
1164   // Initialization (coefficient = alpha, rates = lambda)
1165   Double_t r1 = 1.0;
1166   Double_t r2 = 1.0;
1167   Double_t c1 = 0.5;
1168   Double_t c2 = 0.5;
1169
1170   if (nexp == 1) {   // 1 Exponentials
1171     r1 = 1.156;
1172     r2 = 0.130;
1173     c1 = 0.066;
1174     c2 = 0.000;
1175   }
1176   if (nexp == 2) {   // 2 Exponentials
1177     r1 = 1.156;
1178     r2 = 0.130;
1179     c1 = 0.114;
1180     c2 = 0.624;
1181   }
1182
1183   coefficients[0] = c1;
1184   coefficients[1] = c2;
1185
1186   Double_t dt = 0.1;
1187
1188   rates[0] = TMath::Exp(-dt/(r1));
1189   rates[1] = TMath::Exp(-dt/(r2));
1190
1191   Int_t i = 0;
1192   Int_t k = 0;
1193
1194   Double_t reminder[2];
1195   Double_t correction = 0.0;
1196   Double_t result     = 0.0;
1197
1198   // Attention: computation order is important
1199   for (k = 0; k < nexp; k++) {
1200     reminder[k] = 0.0;
1201   }
1202
1203   for (i = 0; i < n; i++) {
1204
1205     result    = (source[i] - correction);    // No rescaling
1206     target[i] = result;
1207
1208     for (k = 0; k < nexp; k++) {
1209       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1210     }
1211
1212     correction = 0.0;
1213     for (k = 0; k < nexp; k++) {
1214       correction += reminder[k];
1215     }
1216
1217   }
1218
1219 }
1220
1221 //_____________________________________________________________________________
1222 void AliTRDclusterizer::ResetRecPoints()
1223 {
1224   //
1225   // Resets the list of rec points
1226   //
1227
1228   if (fRecPoints) {
1229     fRecPoints->Delete();
1230   }
1231
1232 }
1233
1234 //_____________________________________________________________________________
1235 TObjArray *AliTRDclusterizer::RecPoints()
1236 {
1237   //
1238   // Returns the list of rec points
1239   //
1240
1241   if (!fRecPoints) {
1242     fRecPoints = new TObjArray(400);
1243   }
1244
1245   return fRecPoints;
1246
1247 }
1248
1249 //_____________________________________________________________________________
1250 void AliTRDclusterizer::FillLUT()
1251 {
1252   //
1253   // Create the LUT
1254   //
1255
1256   const Int_t kNlut = 128;
1257
1258   fLUTbin = AliTRDgeometry::kNplan * kNlut;
1259
1260   // The lookup table from Bogdan
1261   Float_t lut[AliTRDgeometry::kNplan][kNlut] = {
1262     {
1263       0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1264       0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1265       0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1266       0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1267       0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1268       0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1269       0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1270       0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1271       0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1272       0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1273       0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1274       0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1275       0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1276       0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1277       0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1278       0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1279     },
1280     {
1281       0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1282       0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1283       0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1284       0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1285       0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1286       0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1287       0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1288       0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1289       0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1290       0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1291       0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1292       0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1293       0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1294       0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1295       0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1296       0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1297     },
1298     {
1299       0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1300       0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1301       0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1302       0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1303       0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1304       0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1305       0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1306       0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1307       0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1308       0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1309       0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1310       0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1311       0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1312       0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1313       0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1314       0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1315     },
1316     {
1317       0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1318       0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1319       0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1320       0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1321       0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1322       0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1323       0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1324       0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1325       0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1326       0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1327       0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1328       0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1329       0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1330       0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1331       0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1332       0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1333     },
1334     {
1335       0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1336       0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1337       0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1338       0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1339       0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1340       0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1341       0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1342       0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1343       0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1344       0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1345       0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1346       0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1347       0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1348       0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1349       0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1350       0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1351     },
1352     {
1353       0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1354       0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1355       0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1356       0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1357       0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1358       0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1359       0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1360       0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1361       0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1362       0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1363       0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1364       0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1365       0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1366       0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1367       0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1368       0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1369     }
1370   };
1371
1372   if (fLUT) {
1373     delete [] fLUT;
1374   }
1375   fLUT = new Double_t[fLUTbin];
1376
1377   for (Int_t iplan = 0; iplan < AliTRDgeometry::kNplan; iplan++) {
1378     for (Int_t ilut  = 0; ilut  <  kNlut; ilut++  ) {
1379       fLUT[iplan*kNlut+ilut] = lut[iplan][ilut];
1380     }
1381   }
1382
1383 }
1384
1385 //_____________________________________________________________________________
1386 Double_t AliTRDclusterizer::LUTposition(Int_t iplane, Double_t ampL
1387                                       , Double_t ampC, Double_t ampR) const
1388 {
1389   //
1390   // Calculates the cluster position using the lookup table.
1391   // Method provided by Bogdan Vulpescu.
1392   //
1393
1394   const Int_t kNlut = 128;
1395
1396   Double_t pos;
1397   Double_t x    = 0.0;
1398   Double_t xmin;
1399   Double_t xmax;
1400   Double_t xwid;
1401
1402   Int_t    side = 0;
1403   Int_t    ix;
1404
1405   Double_t xMin[AliTRDgeometry::kNplan] = { 0.006492, 0.006377, 0.006258
1406                                           , 0.006144, 0.006030, 0.005980 };
1407   Double_t xMax[AliTRDgeometry::kNplan] = { 0.960351, 0.965870, 0.970445
1408                                           , 0.974352, 0.977667, 0.996101 };
1409
1410   if      (ampL > ampR) {
1411     x    = (ampL - ampR) / ampC;
1412     side = -1;
1413   }
1414   else if (ampL < ampR) {
1415     x    = (ampR - ampL) / ampC;
1416     side = +1;
1417   }
1418
1419   if (ampL != ampR) {
1420
1421     xmin = xMin[iplane] + 0.000005;
1422     xmax = xMax[iplane] - 0.000005;
1423     xwid = (xmax - xmin) / 127.0;
1424
1425     if      (x < xmin) {
1426       pos = 0.0000;
1427     }
1428     else if (x > xmax) {
1429       pos = side * 0.5000;
1430     }
1431     else {
1432       ix  = (Int_t) ((x - xmin) / xwid);
1433       pos = side * fLUT[iplane*kNlut+ix];
1434     }
1435
1436   }
1437   else {
1438
1439     pos = 0.0;
1440
1441   }
1442
1443   return pos;
1444
1445 }