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