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