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