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