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