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