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