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