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