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   // Threshold value for the maximum ( cut noise)
 649   Float_t minMaxCutSigma = AliTRDReconstructor::RecoParam()->GetMinMaxCutSigma();
 650   // Threshold value for the sum pad ( cut noise)
 651   Float_t minLeftRightCutSigma = AliTRDReconstructor::RecoParam()->GetMinLeftRightCutSigma();
 652
 653   // Iteration limit for unfolding procedure
 654   const Float_t kEpsilon = 0.01;
 655   const Int_t   kNclus   = 3;
 656   const Int_t   kNsig    = 5;
 657
 658   Int_t    iUnfold       = 0;
 659   Double_t ratioLeft     = 1.0;
 660   Double_t ratioRight    = 1.0;
 661
 662   Double_t padSignal[kNsig];
 663   Double_t clusterSignal[kNclus];
 664
 665   Int_t istack  = indexesIn->GetStack();
 666   Int_t ilayer  = indexesIn->GetLayer();
 667   Int_t isector = indexesIn->GetSM();
 668
 669   // Start clustering in the chamber
 670
 671   Int_t idet  = AliTRDgeometry::GetDetector(ilayer,istack,isector);
 672   if (idet != det)
 673     {
 674       AliError("Strange Detector number Missmatch!");
 675       return kFALSE;
 676     }
 677
 678   // TRD space point transformation
 679   fTransform->SetDetector(det);
 680
 681   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + ilayer;
 682   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 683   UShort_t volid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 684
 685   Int_t nColMax    = digitsIn->GetNcol();
 686   Int_t nRowMax    = digitsIn->GetNrow();
 687   Int_t nTimeTotal = digitsIn->GetNtime();
 688
 689   // Detector wise calibration object for the gain factors
 690   const AliTRDCalDet           *calGainFactorDet      = calibration->GetGainFactorDet();
 691   // Calibration object with pad wise values for the gain factors
 692   AliTRDCalROC                 *calGainFactorROC      = calibration->GetGainFactorROC(idet);
 693   // Calibration value for chamber wise gain factor
 694   Float_t                       calGainFactorDetValue = calGainFactorDet->GetValue(idet);
 695
 696
 697   // Detector wise calibration object for the noise
 698   const AliTRDCalDet           *calNoiseDet           = calibration->GetNoiseDet();
 699   // Calibration object with pad wise values for the noise
 700   AliTRDCalROC                 *calNoiseROC           = calibration->GetNoiseROC(idet);
 701   // Calibration value for chamber wise noise
 702   Float_t                       calNoiseDetValue      = calNoiseDet->GetValue(idet);
 703
 704   Int_t nClusters = 0;
 705
 706   AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(nRowMax, nColMax, nTimeTotal);
 707   AliTRDdataArrayS padStatus(nRowMax, nColMax, nTimeTotal);
 708
 709   ResetHelperIndexes(indexesIn);
 710
 711   // Apply the gain and the tail cancelation via digital filter
 712   TailCancelation(digitsIn
 713                 ,digitsOut
 714                 ,indexesIn
 715                 ,fIndexesOut
 716                 ,nTimeTotal
 717                 ,adcThreshold
 718                 ,calGainFactorROC
 719                 ,calGainFactorDetValue);
 720
 721   Int_t row  = 0;
 722   Int_t col  = 0;
 723   Int_t time = 0;
 724   Int_t iPad = 0;
 725
 726   UChar_t status[3]={0, 0, 0}, ipos = 0;
 727   fIndexesOut->ResetCounters();
 728   while (fIndexesOut->NextRCTbinIndex(row, col, time)){
 729     Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
 730     status[1] = digitsIn->GetPadStatus(row,col,time);
 731     if(status[1]) SETBIT(ipos, AliTRDcluster::kMaskedCenter);
 732
 733     if(signalM < maxThresh) continue;
 734
 735     Float_t  noiseMiddleThresh = minMaxCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
 736     if (signalM < noiseMiddleThresh) continue;
 737
 738     if (col + 1 >= nColMax || col-1 < 0) continue;
 739
 740     Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
 741     status[0] = digitsIn->GetPadStatus(row,col+1,time);
 742     if(status[0]) SETBIT(ipos, AliTRDcluster::kMaskedLeft);
 743
 744     Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
 745     status[2] = digitsIn->GetPadStatus(row,col-1,time);
 746     if(status[2]) SETBIT(ipos, AliTRDcluster::kMaskedRight);
 747
 748     // reject candidates with more than 1 problematic pad
 749     if(ipos == 3 || ipos > 4) continue;
 750
 751     if(!status[1]){ // good central pad
 752       if(!ipos){ // all pads are OK
 753         if ((signalL <= signalM) && (signalR <  signalM)) {
 754           if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
 755             Float_t  noiseSumThresh    = minLeftRightCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
 756             if((signalL+signalR+signalM) >= noiseSumThresh){
 757               // Maximum found, mark the position by a negative signal
 758               digitsOut->SetDataUnchecked(row,col,time,-signalM);
 759               fIndexesMaxima->AddIndexTBin(row,col,time);
 760               padStatus.SetDataUnchecked(row, col, time, ipos);
 761             }
 762           }
 763         }
 764       } else { // one of the neighbouring pads are bad
 765         if(status[0] && signalR < signalM && signalR >= sigThresh){
 766           digitsOut->SetDataUnchecked(row,col,time,-signalM);
 767           digitsOut->SetDataUnchecked(row, col, time+1, 0.);
 768           fIndexesMaxima->AddIndexTBin(row,col,time);
 769           padStatus.SetDataUnchecked(row, col, time, ipos);
 770         } else if(status[2] && signalL <= signalM && signalL >= sigThresh){
 771           digitsOut->SetDataUnchecked(row,col,time,-signalM);
 772           digitsOut->SetDataUnchecked(row, col, time-1, 0.);
 773               fIndexesMaxima->AddIndexTBin(row,col,time);
 774               padStatus.SetDataUnchecked(row, col, time, ipos);
 775         }
 776       }
 777     } else { // wrong maximum pad
 778       if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
 779         // Maximum found, mark the position by a negative signal
 780         digitsOut->SetDataUnchecked(row,col,time,-maxThresh);
 781         fIndexesMaxima->AddIndexTBin(row,col,time);
 782         padStatus.SetDataUnchecked(row, col, time, ipos);
 783       }
 784     }
 785   }
 786
 787     // The index to the first cluster of a given ROC
 788   Int_t firstClusterROC = -1;
 789     // The number of cluster in a given ROC
 790     Int_t nClusterROC     =  0;
 791
 792     // Now check the maxima and calculate the cluster position
 793     fIndexesMaxima->ResetCounters();
 794     while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) {
 795
 796       // Maximum found ?
 797       if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
 798         for (iPad = 0; iPad < kNclus; iPad++) {
 799           Int_t iPadCol = col - 1 + iPad;
 800           clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
 801         }
 802
 803         // Count the number of pads in the cluster
 804         Int_t nPadCount = 0;
 805         Int_t ii;
 806         // Look to the left
 807         ii = 0;
 808         while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii  ,time)) >= sigThresh) {
 809           nPadCount++;
 810           ii++;
 811           if (col-ii   <        0) break;
 812         }
 813         // Look to the right
 814         ii = 0;
 815         while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh){
 816           nPadCount++;
 817           ii++;
 818           if (col+ii+1 >= nColMax) break;
 819         }
 820         nClusters++;
 821
 822         // Look for 5 pad cluster with minimum in the middle
 823         Bool_t fivePadCluster = kFALSE;
 824         if (col < (nColMax - 3)){
 825           if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
 826             fivePadCluster = kTRUE;
 827           }
 828           if ((fivePadCluster) && (col < (nColMax - 5))) {
 829             if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh){
 830               fivePadCluster = kFALSE;
 831             }
 832           }
 833           if ((fivePadCluster) && (col >             1)){
 834             if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh){
 835               fivePadCluster = kFALSE;
 836             }
 837           }
 838         }
 839
 840         // 5 pad cluster
 841         // Modify the signal of the overlapping pad for the left part
 842         // of the cluster which remains from a previous unfolding
 843         if (iUnfold) {
 844           clusterSignal[0] *= ratioLeft;
 845           iUnfold = 0;
 846         }
 847
 848         // Unfold the 5 pad cluster
 849         if (fivePadCluster){
 850           for (iPad = 0; iPad < kNsig; iPad++) {
 851             padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
 852                       ,col-1+iPad
 853                       ,time));
 854           }
 855           // Unfold the two maxima and set the signal on
 856           // the overlapping pad to the ratio
 857           ratioRight        = Unfold(kEpsilon,ilayer,padSignal);
 858           ratioLeft         = 1.0 - ratioRight;
 859           clusterSignal[2] *= ratioRight;
 860           iUnfold = 1;
 861         }
 862
 863         // The position of the cluster in COL direction relative to the center pad (pad units)
 864         Double_t clusterPosCol = 0.0;
 865         if (AliTRDReconstructor::RecoParam()->LUTOn()) {
 866           // Calculate the position of the cluster by using the
 867           // lookup table method
 868           clusterPosCol = LUTposition(ilayer,clusterSignal[0]
 869               ,clusterSignal[1]
 870               ,clusterSignal[2]);
 871         } else {
 872           // Calculate the position of the cluster by using the
 873           // center of gravity method
 874           for (Int_t i = 0; i < kNsig; i++) {
 875             padSignal[i] = 0.0;
 876           }
 877           padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col  ,time)); // Central pad
 878           padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left    pad
 879           padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right   pad
 880
 881           if ((col >           2) &&
 882               (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
 883               padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
 884           }
 885           if ((col < nColMax - 3) &&
 886               (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])){
 887               padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
 888           }
 889           clusterPosCol = GetCOG(padSignal);
 890         }
 891
 892         // Store the amplitudes of the pads in the cluster for later analysis
 893         // and check whether one of these pads is masked in the database
 894         Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
 895         for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
 896           if ((jPad <          0) ||
 897               (jPad >= nColMax-1)) {
 898               continue;
 899           }
 900           signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
 901         }
 902
 903         // Transform the local cluster coordinates into calibrated
 904         // space point positions defined in the local tracking system.
 905         // Here the calibration for T0, Vdrift and ExB is applied as well.
 906         Double_t clusterXYZ[6];
 907         clusterXYZ[0] = clusterPosCol;
 908         clusterXYZ[1] = clusterSignal[0];
 909         clusterXYZ[2] = clusterSignal[1];
 910         clusterXYZ[3] = clusterSignal[2];
 911         clusterXYZ[4] = 0.0;
 912         clusterXYZ[5] = 0.0;
 913         Int_t    clusterRCT[3];
 914         clusterRCT[0] = row;
 915         clusterRCT[1] = col;
 916         clusterRCT[2] = 0;
 917
 918         Bool_t out = kTRUE;
 919         if (fTransform->Transform(clusterXYZ, clusterRCT, ((UInt_t) time), out, 0)) {
 920
 921         // Add the cluster to the output array
 922         // The track indices will be stored later
 923         Float_t clusterPos[3];
 924         clusterPos[0] = clusterXYZ[0];
 925         clusterPos[1] = clusterXYZ[1];
 926         clusterPos[2] = clusterXYZ[2];
 927         Float_t clusterSig[2];
 928         clusterSig[0] = clusterXYZ[4];
 929         clusterSig[1] = clusterXYZ[5];
 930         Double_t clusterCharge  = clusterXYZ[3];
 931         Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
 932         AliTRDcluster *cluster = new AliTRDcluster(idet
 933                     ,clusterCharge
 934                     ,clusterPos
 935                     ,clusterSig
 936                     ,0x0
 937                     ,((Char_t) nPadCount)
 938                     ,signals
 939                     ,((UChar_t) col)
 940                     ,((UChar_t) row)
 941                     ,((UChar_t) time)
 942                     ,clusterTimeBin
 943                     ,clusterPosCol
 944                     ,volid);
 945         cluster->SetInChamber(!out);
 946         UChar_t maskPosition = padStatus.GetDataUnchecked(row, col, time);
 947         if(maskPosition){
 948           cluster->SetPadMaskedPosition(maskPosition);
 949
 950           if(maskPosition & AliTRDcluster::kMaskedLeft) cluster->SetPadMaskedStatus(status[0]);
 951           else if(maskPosition & AliTRDcluster::kMaskedCenter) cluster->SetPadMaskedStatus(status[1]);
 952           else cluster->SetPadMaskedStatus(status[2]);
 953         }
 954
 955         // Temporarily store the row, column and time bin of the center pad
 956         // Used to later on assign the track indices
 957         cluster->SetLabel( row,0);
 958         cluster->SetLabel( col,1);
 959         cluster->SetLabel(time,2);
 960
 961         RecPoints()->Add(cluster);
 962
 963         // Store the index of the first cluster in the current ROC
 964         if (firstClusterROC < 0) {
 965           firstClusterROC = RecPoints()->GetEntriesFast() - 1;
 966         }
 967
 968         // Count the number of cluster in the current ROC
 969         nClusterROC++;
 970
 971       } // if: Transform ok ?
 972     } // if: Maximum found ?
 973   }
 974
 975   delete digitsOut;
 976
 977   if (fAddLabels) AddLabels(idet, firstClusterROC, nClusterROC);
 978
 979   // Write the cluster and reset the array
 980   WriteClusters(idet);
 981   ResetRecPoints();
 982
 983   return kTRUE;
 984
 985 }
 986
 987 //_____________________________________________________________________________
 988 Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
 989 {
 990   //
 991   // Add the track indices to the found clusters
 992   //
 993
 994   const Int_t   kNclus  = 3;
 995   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
 996   const Int_t   kNtrack = kNdict * kNclus;
 997
 998   Int_t iClusterROC = 0;
 999
1000   Int_t row  = 0;
1001   Int_t col  = 0;
1002   Int_t time = 0;
1003   Int_t iPad = 0;
1004
1005   // Temporary array to collect the track indices
1006   Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1007
1008   // Loop through the dictionary arrays one-by-one
1009   // to keep memory consumption low
1010   AliTRDdataArrayI *tracksIn = 0;
1011   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1012
1013     // tracksIn should be expanded beforehand!
1014     tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
1015
1016     // Loop though the clusters found in this ROC
1017     for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1018
1019       AliTRDcluster *cluster = (AliTRDcluster *)
1020         RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1021       row  = cluster->GetLabel(0);
1022       col  = cluster->GetLabel(1);
1023       time = cluster->GetLabel(2);
1024
1025       for (iPad = 0; iPad < kNclus; iPad++) {
1026         Int_t iPadCol = col - 1 + iPad;
1027         Int_t index   = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1028         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1029       }
1030
1031     }
1032
1033   }
1034
1035   // Copy the track indices into the cluster
1036   // Loop though the clusters found in this ROC
1037   for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1038
1039     AliTRDcluster *cluster = (AliTRDcluster *)
1040       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1041     cluster->SetLabel(-9999,0);
1042     cluster->SetLabel(-9999,1);
1043     cluster->SetLabel(-9999,2);
1044
1045     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1046
1047   }
1048
1049   delete [] idxTracks;
1050
1051   return kTRUE;
1052
1053 }
1054
1055 //_____________________________________________________________________________
1056 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
1057 {
1058   //
1059   // Get COG position
1060   // Used for clusters with more than 3 pads - where LUT not applicable
1061   //
1062
1063   Double_t sum = signal[0]
1064                + signal[1]
1065                + signal[2]
1066                + signal[3]
1067                + signal[4];
1068
1069   Double_t res = (0.0 * (-signal[0] + signal[4])
1070                       + (-signal[1] + signal[3])) / sum;
1071
1072   return res;
1073
1074 }
1075
1076 //_____________________________________________________________________________
1077 Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal)
1078 {
1079   //
1080   // Method to unfold neighbouring maxima.
1081   // The charge ratio on the overlapping pad is calculated
1082   // until there is no more change within the range given by eps.
1083   // The resulting ratio is then returned to the calling method.
1084   //
1085
1086   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1087   if (!calibration) {
1088     AliError("No AliTRDcalibDB instance available\n");
1089     return kFALSE;
1090   }
1091
1092   Int_t   irc                = 0;
1093   Int_t   itStep             = 0;                 // Count iteration steps
1094
1095   Double_t ratio             = 0.5;               // Start value for ratio
1096   Double_t prevRatio         = 0.0;               // Store previous ratio
1097
1098   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1099   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1100   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1101
1102   // Start the iteration
1103   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1104
1105     itStep++;
1106     prevRatio = ratio;
1107
1108     // Cluster position according to charge ratio
1109     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1110                       / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1111     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1112                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1113
1114     // Set cluster charge ratio
1115     irc = calibration->PadResponse(1.0,maxLeft ,layer,newSignal);
1116     Double_t ampLeft  = padSignal[1] / newSignal[1];
1117     irc = calibration->PadResponse(1.0,maxRight,layer,newSignal);
1118     Double_t ampRight = padSignal[3] / newSignal[1];
1119
1120     // Apply pad response to parameters
1121     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1122     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1123
1124     // Calculate new overlapping ratio
1125     ratio = TMath::Min((Double_t) 1.0
1126                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1127
1128   }
1129
1130   return ratio;
1131
1132 }
1133
1134 //_____________________________________________________________________________
1135 void AliTRDclusterizer::TailCancelation(AliTRDdataArrayDigits *digitsIn
1136                                       , AliTRDdataArrayF *digitsOut
1137                                       , AliTRDSignalIndex *indexesIn
1138                                       , AliTRDSignalIndex *indexesOut
1139                                       , Int_t nTimeTotal
1140                                       , Float_t adcThreshold
1141                                       , AliTRDCalROC *calGainFactorROC
1142                                       , Float_t calGainFactorDetValue)
1143 {
1144   //
1145   // Applies the tail cancelation and gain factors:
1146   // Transform digitsIn to digitsOut
1147   //
1148
1149   Int_t iRow  = 0;
1150   Int_t iCol  = 0;
1151   Int_t iTime = 0;
1152
1153   Double_t *inADC  = new Double_t[nTimeTotal];  // ADC data before tail cancellation
1154   Double_t *outADC = new Double_t[nTimeTotal];  // ADC data after tail cancellation
1155   indexesIn->ResetCounters();
1156   while (indexesIn->NextRCIndex(iRow, iCol))
1157     {
1158       Float_t  calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1159       Double_t gain                  = calGainFactorDetValue
1160                                      * calGainFactorROCValue;
1161
1162       Bool_t corrupted = kFALSE;
1163       for (iTime = 0; iTime < nTimeTotal; iTime++)
1164         {
1165           // Apply gain gain factor
1166           inADC[iTime]   = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1167           if(digitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
1168           inADC[iTime]  /= gain;
1169           outADC[iTime]  = inADC[iTime];
1170         }
1171       if(!corrupted)
1172         {
1173           // Apply the tail cancelation via the digital filter
1174           // (only for non-coorupted pads)
1175           if (AliTRDReconstructor::RecoParam()->TCOn())
1176             {
1177               DeConvExp(inADC,outADC,nTimeTotal,AliTRDReconstructor::RecoParam()->GetTCnexp());
1178             }
1179         }
1180
1181       indexesIn->ResetTbinCounter();
1182       while (indexesIn->NextTbinIndex(iTime))
1183         {
1184           // Store the amplitude of the digit if above threshold
1185           if (outADC[iTime] > adcThreshold)
1186             {
1187               digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1188               indexesOut->AddIndexTBin(iRow,iCol,iTime);
1189             }
1190         } // while itime
1191
1192     } // while irow icol
1193
1194   delete [] inADC;
1195   delete [] outADC;
1196
1197   return;
1198
1199 }
1200
1201 //_____________________________________________________________________________
1202 void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1203                                  , Int_t n, Int_t nexp)
1204 {
1205   //
1206   // Tail cancellation by deconvolution for PASA v4 TRF
1207   //
1208
1209   Double_t rates[2];
1210   Double_t coefficients[2];
1211
1212   // Initialization (coefficient = alpha, rates = lambda)
1213   Double_t r1 = 1.0;
1214   Double_t r2 = 1.0;
1215   Double_t c1 = 0.5;
1216   Double_t c2 = 0.5;
1217
1218   if (nexp == 1) {   // 1 Exponentials
1219     r1 = 1.156;
1220     r2 = 0.130;
1221     c1 = 0.066;
1222     c2 = 0.000;
1223   }
1224   if (nexp == 2) {   // 2 Exponentials
1225     r1 = 1.156;
1226     r2 = 0.130;
1227     c1 = 0.114;
1228     c2 = 0.624;
1229   }
1230
1231   coefficients[0] = c1;
1232   coefficients[1] = c2;
1233
1234   Double_t dt = 0.1;
1235
1236   rates[0] = TMath::Exp(-dt/(r1));
1237   rates[1] = TMath::Exp(-dt/(r2));
1238
1239   Int_t i = 0;
1240   Int_t k = 0;
1241
1242   Double_t reminder[2];
1243   Double_t correction = 0.0;
1244   Double_t result     = 0.0;
1245
1246   // Attention: computation order is important
1247   for (k = 0; k < nexp; k++) {
1248     reminder[k] = 0.0;
1249   }
1250
1251   for (i = 0; i < n; i++) {
1252
1253     result    = (source[i] - correction);    // No rescaling
1254     target[i] = result;
1255
1256     for (k = 0; k < nexp; k++) {
1257       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1258     }
1259
1260     correction = 0.0;
1261     for (k = 0; k < nexp; k++) {
1262       correction += reminder[k];
1263     }
1264
1265   }
1266
1267 }
1268
1269 //_____________________________________________________________________________
1270 void AliTRDclusterizer::ResetRecPoints()
1271 {
1272   //
1273   // Resets the list of rec points
1274   //
1275
1276   if (fRecPoints) {
1277     fRecPoints->Delete();
1278   }
1279
1280 }
1281
1282 //_____________________________________________________________________________
1283 TObjArray *AliTRDclusterizer::RecPoints()
1284 {
1285   //
1286   // Returns the list of rec points
1287   //
1288
1289   if (!fRecPoints) {
1290     fRecPoints = new TObjArray(400);
1291   }
1292
1293   return fRecPoints;
1294
1295 }
1296
1297 //_____________________________________________________________________________
1298 void AliTRDclusterizer::FillLUT()
1299 {
1300   //
1301   // Create the LUT
1302   //
1303
1304   const Int_t kNlut = 128;
1305
1306   fLUTbin = AliTRDgeometry::kNlayer * kNlut;
1307
1308   // The lookup table from Bogdan
1309   Float_t lut[AliTRDgeometry::kNlayer][kNlut] = {
1310     {
1311       0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1312       0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1313       0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1314       0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1315       0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1316       0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1317       0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1318       0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1319       0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1320       0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1321       0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1322       0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1323       0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1324       0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1325       0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1326       0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1327     },
1328     {
1329       0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1330       0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1331       0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1332       0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1333       0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1334       0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1335       0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1336       0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1337       0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1338       0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1339       0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1340       0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1341       0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1342       0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1343       0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1344       0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1345     },
1346     {
1347       0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1348       0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1349       0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1350       0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1351       0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1352       0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1353       0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1354       0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1355       0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1356       0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1357       0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1358       0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1359       0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1360       0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1361       0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1362       0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1363     },
1364     {
1365       0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1366       0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1367       0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1368       0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1369       0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1370       0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1371       0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1372       0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1373       0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1374       0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1375       0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1376       0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1377       0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1378       0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1379       0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1380       0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1381     },
1382     {
1383       0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1384       0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1385       0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1386       0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1387       0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1388       0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1389       0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1390       0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1391       0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1392       0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1393       0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1394       0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1395       0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1396       0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1397       0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1398       0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1399     },
1400     {
1401       0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1402       0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1403       0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1404       0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1405       0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1406       0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1407       0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1408       0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1409       0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1410       0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1411       0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1412       0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1413       0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1414       0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1415       0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1416       0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1417     }
1418   };
1419
1420   if (fLUT) {
1421     delete [] fLUT;
1422   }
1423   fLUT = new Double_t[fLUTbin];
1424
1425   for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1426     for (Int_t ilut  = 0; ilut  <  kNlut; ilut++  ) {
1427       fLUT[ilayer*kNlut+ilut] = lut[ilayer][ilut];
1428     }
1429   }
1430
1431 }
1432
1433 //_____________________________________________________________________________
1434 Double_t AliTRDclusterizer::LUTposition(Int_t ilayer, Double_t ampL
1435                                       , Double_t ampC, Double_t ampR) const
1436 {
1437   //
1438   // Calculates the cluster position using the lookup table.
1439   // Method provided by Bogdan Vulpescu.
1440   //
1441
1442   const Int_t kNlut = 128;
1443
1444   Double_t pos;
1445   Double_t x    = 0.0;
1446   Double_t xmin;
1447   Double_t xmax;
1448   Double_t xwid;
1449
1450   Int_t    side = 0;
1451   Int_t    ix;
1452
1453   Double_t xMin[AliTRDgeometry::kNlayer] = { 0.006492, 0.006377, 0.006258
1454                                            , 0.006144, 0.006030, 0.005980 };
1455   Double_t xMax[AliTRDgeometry::kNlayer] = { 0.960351, 0.965870, 0.970445
1456                                            , 0.974352, 0.977667, 0.996101 };
1457
1458   if      (ampL > ampR) {
1459     x    = (ampL - ampR) / ampC;
1460     side = -1;
1461   }
1462   else if (ampL < ampR) {
1463     x    = (ampR - ampL) / ampC;
1464     side = +1;
1465   }
1466
1467   if (ampL != ampR) {
1468
1469     xmin = xMin[ilayer] + 0.000005;
1470     xmax = xMax[ilayer] - 0.000005;
1471     xwid = (xmax - xmin) / 127.0;
1472
1473     if      (x < xmin) {
1474       pos = 0.0000;
1475     }
1476     else if (x > xmax) {
1477       pos = side * 0.5000;
1478     }
1479     else {
1480       ix  = (Int_t) ((x - xmin) / xwid);
1481       pos = side * fLUT[ilayer*kNlut+ix];
1482     }
1483
1484   }
1485   else {
1486
1487     pos = 0.0;
1488
1489   }
1490
1491   return pos;
1492
1493 }