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