TRD/AliTRDclusterizer.cxx

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