TRD/AliTRDclusterizer.cxx

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