TRD/AliTRDclusterizer.cxx

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