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