TRD/AliTRDclusterizer.cxx

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