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 <TClonesArray.h>
  25 #include <TObjArray.h>
  26
  27 #include "AliRunLoader.h"
  28 #include "AliLoader.h"
  29 #include "AliAlignObj.h"
  30
  31 #include "AliTRDclusterizer.h"
  32 #include "AliTRDcluster.h"
  33 #include "AliTRDReconstructor.h"
  34 #include "AliTRDgeometry.h"
  35 #include "AliTRDarrayDictionary.h"
  36 #include "AliTRDarrayADC.h"
  37 #include "AliTRDdigitsManager.h"
  38 #include "AliTRDdigitsParam.h"
  39 #include "AliTRDrawData.h"
  40 #include "AliTRDcalibDB.h"
  41 #include "AliTRDtransform.h"
  42 #include "AliTRDSignalIndex.h"
  43 #include "AliTRDrawStreamBase.h"
  44 #include "AliTRDfeeParam.h"
  45 #include "AliTRDtrackletWord.h"
  46
  47 #include "TTreeStream.h"
  48
  49 #include "Cal/AliTRDCalROC.h"
  50 #include "Cal/AliTRDCalDet.h"
  51 #include "Cal/AliTRDCalSingleChamberStatus.h"
  52
  53 ClassImp(AliTRDclusterizer)
  54
  55 //_____________________________________________________________________________
  56 AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
  57   :TNamed()
  58   ,fReconstructor(rec)
  59   ,fRunLoader(NULL)
  60   ,fClusterTree(NULL)
  61   ,fRecPoints(NULL)
  62   ,fTracklets(NULL)
  63   ,fTrackletTree(NULL)
  64   ,fDigitsManager(new AliTRDdigitsManager())
  65   ,fTrackletContainer(NULL)
  66   ,fRawVersion(2)
  67   ,fTransform(new AliTRDtransform(0))
  68   ,fDigits(NULL)
  69   ,fIndexes(NULL)
  70   ,fMaxThresh(0)
  71   ,fSigThresh(0)
  72   ,fMinMaxCutSigma(0)
  73   ,fMinLeftRightCutSigma(0)
  74   ,fLayer(0)
  75   ,fDet(0)
  76   ,fVolid(0)
  77   ,fColMax(0)
  78   ,fTimeTotal(0)
  79   ,fCalGainFactorROC(NULL)
  80   ,fCalGainFactorDetValue(0)
  81   ,fCalNoiseROC(NULL)
  82   ,fCalNoiseDetValue(0)
  83   ,fCalPadStatusROC(NULL)
  84   ,fClusterROC(0)
  85   ,firstClusterROC(0)
  86   ,fNoOfClusters(0)
  87 {
  88   //
  89   // AliTRDclusterizer default constructor
  90   //
  91
  92   SetBit(kLabels, kTRUE);
  93
  94   AliTRDcalibDB *trd = 0x0;
  95   if (!(trd = AliTRDcalibDB::Instance())) {
  96     AliFatal("Could not get calibration object");
  97   }
  98
  99   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 100
 101   // Initialize debug stream
 102   if(fReconstructor){
 103     if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 1){
 104       TDirectory *savedir = gDirectory;
 105       //fgGetDebugStream    = new TTreeSRedirector("TRD.ClusterizerDebug.root");
 106       savedir->cd();
 107     }
 108   }
 109
 110 }
 111
 112 //_____________________________________________________________________________
 113 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, const AliTRDReconstructor *const rec)
 114   :TNamed(name,title)
 115   ,fReconstructor(rec)
 116   ,fRunLoader(NULL)
 117   ,fClusterTree(NULL)
 118   ,fRecPoints(NULL)
 119   ,fTracklets(NULL)
 120   ,fTrackletTree(NULL)
 121   ,fDigitsManager(new AliTRDdigitsManager())
 122   ,fTrackletContainer(NULL)
 123   ,fRawVersion(2)
 124   ,fTransform(new AliTRDtransform(0))
 125   ,fDigits(NULL)
 126   ,fIndexes(NULL)
 127   ,fMaxThresh(0)
 128   ,fSigThresh(0)
 129   ,fMinMaxCutSigma(0)
 130   ,fMinLeftRightCutSigma(0)
 131   ,fLayer(0)
 132   ,fDet(0)
 133   ,fVolid(0)
 134   ,fColMax(0)
 135   ,fTimeTotal(0)
 136   ,fCalGainFactorROC(NULL)
 137   ,fCalGainFactorDetValue(0)
 138   ,fCalNoiseROC(NULL)
 139   ,fCalNoiseDetValue(0)
 140   ,fCalPadStatusROC(NULL)
 141   ,fClusterROC(0)
 142   ,firstClusterROC(0)
 143   ,fNoOfClusters(0)
 144 {
 145   //
 146   // AliTRDclusterizer constructor
 147   //
 148
 149   SetBit(kLabels, kTRUE);
 150
 151   AliTRDcalibDB *trd = 0x0;
 152   if (!(trd = AliTRDcalibDB::Instance())) {
 153     AliFatal("Could not get calibration object");
 154   }
 155
 156   fDigitsManager->CreateArrays();
 157
 158   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 159
 160   //FillLUT();
 161
 162 }
 163
 164 //_____________________________________________________________________________
 165 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
 166   :TNamed(c)
 167   ,fReconstructor(c.fReconstructor)
 168   ,fRunLoader(NULL)
 169   ,fClusterTree(NULL)
 170   ,fRecPoints(NULL)
 171   ,fTracklets(NULL)
 172   ,fTrackletTree(NULL)
 173   ,fDigitsManager(NULL)
 174   ,fTrackletContainer(NULL)
 175   ,fRawVersion(2)
 176   ,fTransform(NULL)
 177   ,fDigits(NULL)
 178   ,fIndexes(NULL)
 179   ,fMaxThresh(0)
 180   ,fSigThresh(0)
 181   ,fMinMaxCutSigma(0)
 182   ,fMinLeftRightCutSigma(0)
 183   ,fLayer(0)
 184   ,fDet(0)
 185   ,fVolid(0)
 186   ,fColMax(0)
 187   ,fTimeTotal(0)
 188   ,fCalGainFactorROC(NULL)
 189   ,fCalGainFactorDetValue(0)
 190   ,fCalNoiseROC(NULL)
 191   ,fCalNoiseDetValue(0)
 192   ,fCalPadStatusROC(NULL)
 193   ,fClusterROC(0)
 194   ,firstClusterROC(0)
 195   ,fNoOfClusters(0)
 196 {
 197   //
 198   // AliTRDclusterizer copy constructor
 199   //
 200
 201   SetBit(kLabels, kTRUE);
 202
 203   //FillLUT();
 204
 205 }
 206
 207 //_____________________________________________________________________________
 208 AliTRDclusterizer::~AliTRDclusterizer()
 209 {
 210   //
 211   // AliTRDclusterizer destructor
 212   //
 213
 214   if (fRecPoints/* && IsClustersOwner()*/){
 215     fRecPoints->Delete();
 216     delete fRecPoints;
 217   }
 218
 219   if (fTracklets){
 220     fTracklets->Delete();
 221     delete fTracklets;
 222   }
 223
 224   if (fDigitsManager) {
 225     delete fDigitsManager;
 226     fDigitsManager = NULL;
 227   }
 228
 229   if (fTrackletContainer){
 230     delete fTrackletContainer;
 231     fTrackletContainer = NULL;
 232   }
 233
 234   if (fTransform){
 235     delete fTransform;
 236     fTransform     = NULL;
 237   }
 238
 239 }
 240
 241 //_____________________________________________________________________________
 242 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
 243 {
 244   //
 245   // Assignment operator
 246   //
 247
 248   if (this != &c)
 249     {
 250       ((AliTRDclusterizer &) c).Copy(*this);
 251     }
 252
 253   return *this;
 254
 255 }
 256
 257 //_____________________________________________________________________________
 258 void AliTRDclusterizer::Copy(TObject &c) const
 259 {
 260   //
 261   // Copy function
 262   //
 263
 264   ((AliTRDclusterizer &) c).fClusterTree   = NULL;
 265   ((AliTRDclusterizer &) c).fRecPoints     = NULL;
 266   ((AliTRDclusterizer &) c).fTrackletTree  = NULL;
 267   ((AliTRDclusterizer &) c).fDigitsManager = NULL;
 268   ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
 269   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
 270   ((AliTRDclusterizer &) c).fTransform     = NULL;
 271   ((AliTRDclusterizer &) c).fDigits      = NULL;
 272   ((AliTRDclusterizer &) c).fIndexes       = NULL;
 273   ((AliTRDclusterizer &) c).fMaxThresh     = 0;
 274   ((AliTRDclusterizer &) c).fSigThresh     = 0;
 275   ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0;
 276   ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0;
 277   ((AliTRDclusterizer &) c).fLayer         = 0;
 278   ((AliTRDclusterizer &) c).fDet           = 0;
 279   ((AliTRDclusterizer &) c).fVolid         = 0;
 280   ((AliTRDclusterizer &) c).fColMax        = 0;
 281   ((AliTRDclusterizer &) c).fTimeTotal     = 0;
 282   ((AliTRDclusterizer &) c).fCalGainFactorROC = NULL;
 283   ((AliTRDclusterizer &) c).fCalGainFactorDetValue = 0;
 284   ((AliTRDclusterizer &) c).fCalNoiseROC   = NULL;
 285   ((AliTRDclusterizer &) c).fCalNoiseDetValue = 0;
 286   ((AliTRDclusterizer &) c).fCalPadStatusROC = NULL;
 287   ((AliTRDclusterizer &) c).fClusterROC    = 0;
 288   ((AliTRDclusterizer &) c).firstClusterROC= 0;
 289   ((AliTRDclusterizer &) c).fNoOfClusters  = 0;
 290 }
 291
 292 //_____________________________________________________________________________
 293 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
 294 {
 295   //
 296   // Opens the AliROOT file. Output and input are in the same file
 297   //
 298
 299   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 300   fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 301
 302   if (!fRunLoader) {
 303     fRunLoader = AliRunLoader::Open(name);
 304   }
 305
 306   if (!fRunLoader) {
 307     AliError(Form("Can not open session for file %s.",name));
 308     return kFALSE;
 309   }
 310
 311   OpenInput(nEvent);
 312   OpenOutput();
 313
 314   return kTRUE;
 315
 316 }
 317
 318 //_____________________________________________________________________________
 319 Bool_t AliTRDclusterizer::OpenOutput()
 320 {
 321   //
 322   // Open the output file
 323   //
 324
 325   if (!fReconstructor->IsWritingClusters()) return kTRUE;
 326
 327   TObjArray *ioArray = 0x0;
 328
 329   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 330   loader->MakeTree("R");
 331
 332   fClusterTree = loader->TreeR();
 333   fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 334
 335   return kTRUE;
 336
 337 }
 338
 339 //_____________________________________________________________________________
 340 Bool_t AliTRDclusterizer::OpenOutput(TTree *const clusterTree)
 341 {
 342   //
 343   // Connect the output tree
 344   //
 345
 346   // clusters writing
 347   if (fReconstructor->IsWritingClusters()){
 348     TObjArray *ioArray = 0x0;
 349     fClusterTree = clusterTree;
 350     fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 351   }
 352   return kTRUE;
 353 }
 354
 355 //_____________________________________________________________________________
 356 Bool_t AliTRDclusterizer::OpenTrackletOutput()
 357 {
 358   //
 359   // Tracklet writing
 360   //
 361
 362   if (fReconstructor->IsWritingTracklets()){
 363     TString evfoldname = AliConfig::GetDefaultEventFolderName();
 364     fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 365
 366     if (!fRunLoader) {
 367       fRunLoader = AliRunLoader::Open("galice.root");
 368     }
 369     if (!fRunLoader) {
 370       AliError(Form("Can not open session for file galice.root."));
 371       return kFALSE;
 372     }
 373
 374     UInt_t **leaves = new UInt_t *[2];
 375     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 376     if (!dl) {
 377       AliError("Could not get the tracklets data loader!");
 378       dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
 379       fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
 380     }
 381     fTrackletTree = dl->Tree();
 382     if (!fTrackletTree)
 383       {
 384        dl->MakeTree();
 385        fTrackletTree = dl->Tree();
 386       }
 387     TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 388     if (!trkbranch)
 389       fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 390   }
 391
 392   return kTRUE;
 393 }
 394
 395 //_____________________________________________________________________________
 396 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
 397 {
 398   //
 399   // Opens a ROOT-file with TRD-hits and reads in the digits-tree
 400   //
 401
 402   // Import the Trees for the event nEvent in the file
 403   fRunLoader->GetEvent(nEvent);
 404
 405   return kTRUE;
 406
 407 }
 408
 409 //_____________________________________________________________________________
 410 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
 411 {
 412   //
 413   // Fills TRDcluster branch in the tree with the clusters
 414   // found in detector = det. For det=-1 writes the tree.
 415   //
 416
 417   if ((det <                      -1) ||
 418       (det >= AliTRDgeometry::Ndet())) {
 419     AliError(Form("Unexpected detector index %d.\n",det));
 420     return kFALSE;
 421   }
 422
 423   TObjArray *ioArray = new TObjArray(400);
 424   TBranch *branch = fClusterTree->GetBranch("TRDcluster");
 425   if (!branch) {
 426     branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
 427   } else branch->SetAddress(&ioArray);
 428
 429   Int_t nRecPoints = RecPoints()->GetEntriesFast();
 430   if(det >= 0){
 431     for (Int_t i = 0; i < nRecPoints; i++) {
 432       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 433       if(det != c->GetDetector()) continue;
 434       ioArray->AddLast(c);
 435     }
 436     fClusterTree->Fill();
 437   } else {
 438
 439     Int_t detOld = -1;
 440     for (Int_t i = 0; i < nRecPoints; i++) {
 441       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 442       if(c->GetDetector() != detOld){
 443         fClusterTree->Fill();
 444         ioArray->Clear();
 445         detOld = c->GetDetector();
 446       }
 447       ioArray->AddLast(c);
 448     }
 449   }
 450   delete ioArray;
 451
 452   return kTRUE;
 453
 454 }
 455
 456 //_____________________________________________________________________________
 457 Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
 458 {
 459   //
 460   // Write the raw data tracklets into seperate file
 461   //
 462
 463   UInt_t **leaves = new UInt_t *[2];
 464   for (Int_t i=0; i<2 ;i++){
 465     leaves[i] = new UInt_t[258];
 466     leaves[i][0] = det; // det
 467     leaves[i][1] = i;   // side
 468     memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
 469   }
 470
 471   if (!fTrackletTree){
 472     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 473     dl->MakeTree();
 474     fTrackletTree = dl->Tree();
 475   }
 476
 477   TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 478   if (!trkbranch) {
 479     trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 480   }
 481
 482   for (Int_t i=0; i<2; i++){
 483     if (leaves[i][2]>0) {
 484       trkbranch->SetAddress(leaves[i]);
 485       fTrackletTree->Fill();
 486     }
 487   }
 488
 489   AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 490   dl->WriteData("OVERWRITE");
 491   //dl->Unload();
 492   delete [] leaves;
 493
 494   return kTRUE;
 495
 496 }
 497
 498 //_____________________________________________________________________________
 499 Bool_t AliTRDclusterizer::ReadDigits()
 500 {
 501   //
 502   // Reads the digits arrays from the input aliroot file
 503   //
 504
 505   if (!fRunLoader) {
 506     AliError("No run loader available");
 507     return kFALSE;
 508   }
 509
 510   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 511   if (!loader->TreeD()) {
 512     loader->LoadDigits();
 513   }
 514
 515   // Read in the digit arrays
 516   return (fDigitsManager->ReadDigits(loader->TreeD()));
 517
 518 }
 519
 520 //_____________________________________________________________________________
 521 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
 522 {
 523   //
 524   // Reads the digits arrays from the input tree
 525   //
 526
 527   // Read in the digit arrays
 528   return (fDigitsManager->ReadDigits(digitsTree));
 529
 530 }
 531
 532 //_____________________________________________________________________________
 533 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
 534 {
 535   //
 536   // Reads the digits arrays from the ddl file
 537   //
 538
 539   AliTRDrawData raw;
 540   fDigitsManager = raw.Raw2Digits(rawReader);
 541
 542   return kTRUE;
 543
 544 }
 545
 546 //_____________________________________________________________________________
 547 Bool_t AliTRDclusterizer::MakeClusters()
 548 {
 549   //
 550   // Creates clusters from digits
 551   //
 552
 553   // Propagate info from the digits manager
 554   if (TestBit(kLabels)){
 555     SetBit(kLabels, fDigitsManager->UsesDictionaries());
 556   }
 557
 558   Bool_t fReturn = kTRUE;
 559   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
 560
 561     AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
 562     // This is to take care of switched off super modules
 563     if (!digitsIn->HasData()) continue;
 564     digitsIn->Expand();
 565     digitsIn->DeleteNegatives();  // Restore digits array to >=0 values
 566     AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 567     if (indexes->IsAllocated() == kFALSE){
 568       fDigitsManager->BuildIndexes(i);
 569     }
 570
 571     Bool_t fR = kFALSE;
 572     if (indexes->HasEntry()){
 573       if (TestBit(kLabels)){
 574         for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
 575           AliTRDarrayDictionary *tracksIn = 0; //mod
 576           tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict);  //mod
 577           tracksIn->Expand();
 578         }
 579       }
 580       fR = MakeClusters(i);
 581       fReturn = fR && fReturn;
 582     }
 583
 584     //if (fR == kFALSE){
 585     //  if(IsWritingClusters()) WriteClusters(i);
 586     //  ResetRecPoints();
 587     //}
 588
 589     // No compress just remove
 590     fDigitsManager->RemoveDigits(i);
 591     fDigitsManager->RemoveDictionaries(i);
 592     fDigitsManager->ClearIndexes(i);
 593   }
 594
 595   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 596
 597   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
 598
 599   return fReturn;
 600
 601 }
 602
 603 //_____________________________________________________________________________
 604 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
 605 {
 606   //
 607   // Creates clusters from raw data
 608   //
 609
 610   return Raw2ClustersChamber(rawReader);
 611
 612 }
 613
 614 //_____________________________________________________________________________
 615 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 616 {
 617   //
 618   // Creates clusters from raw data
 619   //
 620
 621   // Create the digits manager
 622   if (!fDigitsManager){
 623     fDigitsManager = new AliTRDdigitsManager(kTRUE);
 624     fDigitsManager->CreateArrays();
 625   }
 626
 627   fDigitsManager->SetUseDictionaries(TestBit(kLabels));
 628
 629   // tracklet container for raw tracklet writing
 630   if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
 631     // maximum tracklets for one HC
 632     const Int_t kTrackletChmb=256;
 633     fTrackletContainer = new UInt_t *[2];
 634     fTrackletContainer[0] = new UInt_t[kTrackletChmb];
 635     fTrackletContainer[1] = new UInt_t[kTrackletChmb];
 636   }
 637
 638   AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader);
 639   if(fReconstructor->IsHLT())
 640     input->SetSharedPadReadout(kFALSE);
 641
 642   AliInfo(Form("Stream version: %s", input->IsA()->GetName()));
 643
 644   Int_t det    = 0;
 645   while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
 646     Bool_t iclusterBranch = kFALSE;
 647     if (fDigitsManager->GetIndexes(det)->HasEntry()){
 648       iclusterBranch = MakeClusters(det);
 649     }
 650
 651     fDigitsManager->ClearArrays(det);
 652
 653     if (!fReconstructor->IsWritingTracklets()) continue;
 654     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
 655   }
 656
 657   if (fTrackletContainer){
 658     delete [] fTrackletContainer[0];
 659     delete [] fTrackletContainer[1];
 660     delete [] fTrackletContainer;
 661     fTrackletContainer = NULL;
 662   }
 663
 664   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 665
 666   delete fDigitsManager;
 667   fDigitsManager = NULL;
 668
 669   delete input;
 670   input = NULL;
 671
 672   AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
 673   return kTRUE;
 674
 675 }
 676
 677 //_____________________________________________________________________________
 678 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 679 {
 680   //
 681   // Check if a pad is masked
 682   //
 683
 684   UChar_t status = 0;
 685
 686   if(signal>0 && TESTBIT(signal, 10)){
 687     CLRBIT(signal, 10);
 688     for(int ibit=0; ibit<4; ibit++){
 689       if(TESTBIT(signal, 11+ibit)){
 690         SETBIT(status, ibit);
 691         CLRBIT(signal, 11+ibit);
 692       }
 693     }
 694   }
 695   return status;
 696 }
 697
 698 //_____________________________________________________________________________
 699 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
 700   //
 701   // Set the pad status into out
 702   // First three bits are needed for the position encoding
 703   //
 704   out |= status << 3;
 705 }
 706
 707 //_____________________________________________________________________________
 708 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
 709   //
 710   // return the staus encoding of the corrupted pad
 711   //
 712   return static_cast<UChar_t>(encoding >> 3);
 713 }
 714
 715 //_____________________________________________________________________________
 716 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
 717   //
 718   // Return the position of the corruption
 719   //
 720   return encoding & 7;
 721 }
 722
 723 //_____________________________________________________________________________
 724 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 725 {
 726   //
 727   // Generates the cluster.
 728   //
 729
 730   // Get the digits
 731   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
 732
 733   // This is to take care of switched off super modules
 734   if (!fDigits->HasData()) return kFALSE;
 735
 736   fIndexes = fDigitsManager->GetIndexes(det);
 737   if (fIndexes->IsAllocated() == kFALSE) {
 738     AliError("Indexes do not exist!");
 739     return kFALSE;
 740   }
 741
 742   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
 743   if (!calibration) {
 744     AliFatal("No AliTRDcalibDB instance available\n");
 745     return kFALSE;
 746   }
 747
 748   if (!fReconstructor){
 749     AliError("Reconstructor not set\n");
 750     return kFALSE;
 751   }
 752
 753   fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
 754   fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
 755   fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
 756   fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
 757
 758   Int_t istack  = fIndexes->GetStack();
 759   fLayer  = fIndexes->GetLayer();
 760   Int_t isector = fIndexes->GetSM();
 761
 762   // Start clustering in the chamber
 763
 764   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 765   if (fDet != det) {
 766     AliError("Strange Detector number Missmatch!");
 767     return kFALSE;
 768   }
 769
 770   // TRD space point transformation
 771   fTransform->SetDetector(det);
 772
 773   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 774   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 775   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 776
 777   if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
 778     AddTrackletsToArray();
 779
 780   fColMax    = fDigits->GetNcol();
 781   //Int_t nRowMax    = fDigits->GetNrow();
 782   fTimeTotal = fDigits->GetNtime();
 783
 784   // Detector wise calibration object for the gain factors
 785   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 786   // Calibration object with pad wise values for the gain factors
 787   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 788   // Calibration value for chamber wise gain factor
 789   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 790
 791   // Detector wise calibration object for the noise
 792   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 793   // Calibration object with pad wise values for the noise
 794   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 795   // Calibration value for chamber wise noise
 796   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 797
 798   // Calibration object with the pad status
 799   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
 800
 801   SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
 802   SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
 803   SetBit(kHLT, fReconstructor->IsHLT());
 804
 805   firstClusterROC = -1;
 806   fClusterROC     =  0;
 807
 808   // Apply the gain and the tail cancelation via digital filter
 809   if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
 810
 811   MaxStruct curr, last;
 812   Int_t nMaximas = 0, nCorrupted = 0;
 813
 814   // Here the clusterfining is happening
 815
 816   for(curr.Time = 0; curr.Time < fTimeTotal; curr.Time++){
 817     while(fIndexes->NextRCIndex(curr.Row, curr.Col)){
 818       //printf("\nCHECK r[%2d] c[%3d] t[%d]\n", curr.Row, curr.Col, curr.Time);
 819       if(IsMaximum(curr, curr.padStatus, &curr.Signals[0])){
 820         //printf("\tMAX s[%d %d %d]\n", curr.Signals[0], curr.Signals[1], curr.Signals[2]);
 821         if(last.Row>-1){
 822           if(curr.Time==last.Time && curr.Row==last.Row && curr.Col==last.Col+2) FivePadCluster(last, curr);
 823           CreateCluster(last);
 824         }
 825         last=curr; curr.FivePad=kFALSE;
 826       }
 827       //printf("\t--- s[%d %d %d]\n", curr.Signals[0], curr.Signals[1], curr.Signals[2]);
 828     }
 829   }
 830   if(last.Row>-1) CreateCluster(last);
 831
 832   if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
 833     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
 834     (*fDebugStream) << "MakeClusters"
 835       << "Detector="   << det
 836       << "NMaxima="    << nMaximas
 837       << "NClusters="  << fClusterROC
 838       << "NCorrupted=" << nCorrupted
 839       << "\n";
 840   }
 841   if (TestBit(kLabels)) AddLabels();
 842
 843   return kTRUE;
 844
 845 }
 846
 847 //_____________________________________________________________________________
 848 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
 849 {
 850   //
 851   // Returns true if this row,col,time combination is a maximum.
 852   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 853   //
 854
 855   Signals[1] = fDigits->GetData(Max.Row, Max.Col, Max.Time);
 856   if(Signals[1] < fMaxThresh) return kFALSE;
 857
 858   Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.Col, Max.Row);
 859   if (Signals[1] < noiseMiddleThresh) return kFALSE;
 860
 861   if (Max.Col + 1 >= fColMax || Max.Col < 1) return kFALSE;
 862
 863   UChar_t status[3]={
 864     fCalPadStatusROC->GetStatus(Max.Col-1, Max.Row)
 865    ,fCalPadStatusROC->GetStatus(Max.Col,   Max.Row)
 866    ,fCalPadStatusROC->GetStatus(Max.Col+1, Max.Row)
 867   };
 868
 869   Signals[0] = fDigits->GetData(Max.Row, Max.Col-1, Max.Time);
 870   Signals[2] = fDigits->GetData(Max.Row, Max.Col+1, Max.Time);
 871
 872   if(!(status[0] | status[1] | status[2])) {//all pads are good
 873     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 874       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
 875         Float_t  noiseSumThresh = fMinLeftRightCutSigma
 876           * fCalNoiseDetValue
 877           * fCalNoiseROC->GetValue(Max.Col, Max.Row);
 878         if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
 879         padStatus = 0;
 880         return kTRUE;
 881       }
 882     }
 883   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 884     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
 885       Signals[2]=0;
 886       SetPadStatus(status[2], padStatus);
 887       return kTRUE;
 888     }
 889     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
 890       Signals[0]=0;
 891       SetPadStatus(status[0], padStatus);
 892       return kTRUE;
 893     }
 894     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
 895       Signals[1]=TMath::Nint(fMaxThresh);
 896       SetPadStatus(status[1], padStatus);
 897       return kTRUE;
 898     }
 899   }
 900   return kFALSE;
 901 }
 902
 903 //_____________________________________________________________________________
 904 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
 905 {
 906   //
 907   // Look for 5 pad cluster with minimum in the middle
 908   // Gives back the ratio
 909   //
 910   if (ThisMax.Col >= fColMax - 3) return kFALSE;
 911   if (ThisMax.Col < fColMax - 5){
 912     if (fDigits->GetData(ThisMax.Row, ThisMax.Col+4, ThisMax.Time) >= fSigThresh)
 913       return kFALSE;
 914   }
 915   if (ThisMax.Col > 1) {
 916     if (fDigits->GetData(ThisMax.Row, ThisMax.Col-2, ThisMax.Time) >= fSigThresh)
 917       return kFALSE;
 918   }
 919
 920   const Float_t kEpsilon = 0.01;
 921   Double_t padSignal[5] = {ThisMax.Signals[0], ThisMax.Signals[1], ThisMax.Signals[2],
 922       NeighbourMax.Signals[1], NeighbourMax.Signals[2]};
 923
 924   // Unfold the two maxima and set the signal on
 925   // the overlapping pad to the ratio
 926   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
 927   ThisMax.Signals[2] = TMath::Nint(ThisMax.Signals[2]*ratio);
 928   NeighbourMax.Signals[0] = TMath::Nint(NeighbourMax.Signals[0]*(1-ratio));
 929   ThisMax.FivePad=kTRUE;
 930   NeighbourMax.FivePad=kTRUE;
 931   return kTRUE;
 932
 933 }
 934
 935 //_____________________________________________________________________________
 936 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 937 {
 938   //
 939   // Creates a cluster at the given position and saves it in fRecPoints
 940   //
 941
 942   Int_t nPadCount = 1;
 943   Short_t signals[7] = { 0, 0, Max.Signals[0], Max.Signals[1], Max.Signals[2], 0, 0 };
 944   if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
 945
 946   AliTRDcluster cluster(fDet, ((UChar_t) Max.Col), ((UChar_t) Max.Row), ((UChar_t) Max.Time), signals, fVolid);
 947   cluster.SetNPads(nPadCount);
 948   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
 949   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
 950   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
 951
 952   cluster.SetFivePad(Max.FivePad);
 953   // set pads status for the cluster
 954   UChar_t maskPosition = GetCorruption(Max.padStatus);
 955   if (maskPosition) {
 956     cluster.SetPadMaskedPosition(maskPosition);
 957     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
 958   }
 959
 960   // Transform the local cluster coordinates into calibrated
 961   // space point positions defined in the local tracking system.
 962   // Here the calibration for T0, Vdrift and ExB is applied as well.
 963   if(!fTransform->Transform(&cluster)) return;
 964   // Temporarily store the Max.Row, column and time bin of the center pad
 965   // Used to later on assign the track indices
 966   cluster.SetLabel(Max.Row, 0);
 967   cluster.SetLabel(Max.Col, 1);
 968   cluster.SetLabel(Max.Time,2);
 969
 970   //needed for HLT reconstruction
 971   AddClusterToArray(&cluster);
 972
 973   // Store the index of the first cluster in the current ROC
 974   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
 975
 976   fNoOfClusters++;
 977   fClusterROC++;
 978 }
 979
 980 //_____________________________________________________________________________
 981 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
 982 {
 983   // Look to the right
 984   Int_t ii = 1;
 985   while (fDigits->GetData(Max.Row, Max.Col-ii, Max.Time) >= fSigThresh) {
 986     nPadCount++;
 987     ii++;
 988     if (Max.Col < ii) break;
 989   }
 990   // Look to the left
 991   ii = 1;
 992   while (fDigits->GetData(Max.Row, Max.Col+ii, Max.Time) >= fSigThresh) {
 993     nPadCount++;
 994     ii++;
 995     if (Max.Col+ii >= fColMax) break;
 996   }
 997
 998   // Store the amplitudes of the pads in the cluster for later analysis
 999   // and check whether one of these pads is masked in the database
1000   signals[2]=Max.Signals[0];
1001   signals[3]=Max.Signals[1];
1002   signals[4]=Max.Signals[2];
1003   for(Int_t i = 0; i<2; i++)
1004     {
1005       if(Max.Col+i >= 3)
1006         signals[i] = fDigits->GetData(Max.Row, Max.Col-3+i, Max.Time);
1007       if(Max.Col+3-i < fColMax)
1008         signals[6-i] = fDigits->GetData(Max.Row, Max.Col+3-i, Max.Time);
1009     }
1010   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1011     if ((jPad >= 0) && (jPad < fColMax))
1012       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1013       }*/
1014 }
1015
1016 //_____________________________________________________________________________
1017 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1018 {
1019   //
1020   // Add a cluster to the array
1021   //
1022
1023   Int_t n = RecPoints()->GetEntriesFast();
1024   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1025   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1026 }
1027
1028 //_____________________________________________________________________________
1029 void AliTRDclusterizer::AddTrackletsToArray()
1030 {
1031   //
1032   // Add the online tracklets of this chamber to the array
1033   //
1034
1035   UInt_t* trackletword;
1036   for(Int_t side=0; side<2; side++)
1037     {
1038       Int_t trkl=0;
1039       trackletword=fTrackletContainer[side];
1040       while(trackletword[trkl]>0){
1041         Int_t n = TrackletsArray()->GetEntriesFast();
1042         AliTRDtrackletWord tmp(trackletword[trkl]);
1043         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1044         trkl++;
1045       }
1046     }
1047 }
1048
1049 //_____________________________________________________________________________
1050 Bool_t AliTRDclusterizer::AddLabels()
1051 {
1052   //
1053   // Add the track indices to the found clusters
1054   //
1055
1056   const Int_t   kNclus  = 3;
1057   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1058   const Int_t   kNtrack = kNdict * kNclus;
1059
1060   Int_t iClusterROC = 0;
1061
1062   Int_t row  = 0;
1063   Int_t col  = 0;
1064   Int_t time = 0;
1065   Int_t iPad = 0;
1066
1067   // Temporary array to collect the track indices
1068   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1069
1070   // Loop through the dictionary arrays one-by-one
1071   // to keep memory consumption low
1072   AliTRDarrayDictionary *tracksIn = 0;  //mod
1073   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1074
1075     // tracksIn should be expanded beforehand!
1076     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1077
1078     // Loop though the clusters found in this ROC
1079     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1080
1081       AliTRDcluster *cluster = (AliTRDcluster *)
1082                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1083       row  = cluster->GetLabel(0);
1084       col  = cluster->GetLabel(1);
1085       time = cluster->GetLabel(2);
1086
1087       for (iPad = 0; iPad < kNclus; iPad++) {
1088         Int_t iPadCol = col - 1 + iPad;
1089         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1090         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1091       }
1092
1093     }
1094
1095   }
1096
1097   // Copy the track indices into the cluster
1098   // Loop though the clusters found in this ROC
1099   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1100
1101     AliTRDcluster *cluster = (AliTRDcluster *)
1102       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1103     cluster->SetLabel(-9999,0);
1104     cluster->SetLabel(-9999,1);
1105     cluster->SetLabel(-9999,2);
1106
1107     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1108
1109   }
1110
1111   delete [] idxTracks;
1112
1113   return kTRUE;
1114
1115 }
1116
1117 //_____________________________________________________________________________
1118 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1119 {
1120   //
1121   // Method to unfold neighbouring maxima.
1122   // The charge ratio on the overlapping pad is calculated
1123   // until there is no more change within the range given by eps.
1124   // The resulting ratio is then returned to the calling method.
1125   //
1126
1127   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1128   if (!calibration) {
1129     AliError("No AliTRDcalibDB instance available\n");
1130     return kFALSE;
1131   }
1132
1133   Int_t   irc                = 0;
1134   Int_t   itStep             = 0;                 // Count iteration steps
1135
1136   Double_t ratio             = 0.5;               // Start value for ratio
1137   Double_t prevRatio         = 0.0;               // Store previous ratio
1138
1139   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1140   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1141   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1142
1143   // Start the iteration
1144   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1145
1146     itStep++;
1147     prevRatio = ratio;
1148
1149     // Cluster position according to charge ratio
1150     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1151                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1152     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1153                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1154
1155     // Set cluster charge ratio
1156     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1157     Double_t ampLeft  = padSignal[1] / newSignal[1];
1158     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1159     Double_t ampRight = padSignal[3] / newSignal[1];
1160
1161     // Apply pad response to parameters
1162     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1163     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1164
1165     // Calculate new overlapping ratio
1166     ratio = TMath::Min((Double_t) 1.0
1167                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1168
1169   }
1170
1171   return ratio;
1172
1173 }
1174
1175 //_____________________________________________________________________________
1176 void AliTRDclusterizer::TailCancelation()
1177 {
1178   //
1179   // Applies the tail cancelation and gain factors:
1180   // Transform fDigits to fDigits
1181   //
1182
1183   Int_t iRow  = 0;
1184   Int_t iCol  = 0;
1185   Int_t iTime = 0;
1186
1187   Double_t *inADC = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1188   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1189
1190   fIndexes->ResetCounters();
1191   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1192   while(fIndexes->NextRCIndex(iRow, iCol))
1193     {
1194       Float_t  fCalGainFactorROCValue = fCalGainFactorROC->GetValue(iCol,iRow);
1195       Double_t gain                  = fCalGainFactorDetValue
1196                                     * fCalGainFactorROCValue;
1197
1198       Bool_t corrupted = kFALSE;
1199       for (iTime = 0; iTime < fTimeTotal; iTime++)
1200         {
1201           // Apply gain gain factor
1202           inADC[iTime]   = fDigits->GetData(iRow,iCol,iTime);
1203           if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1204           inADC[iTime]  /= gain;
1205           outADC[iTime]  = inADC[iTime];
1206           if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming()){
1207             (*fDebugStream) << "TailCancellation"
1208                               << "col="  << iCol
1209                               << "row="  << iRow
1210                               << "time=" << iTime
1211                               << "inADC=" << inADC[iTime]
1212                               << "gain=" << gain
1213                               << "outADC=" << outADC[iTime]
1214                               << "corrupted=" << corrupted
1215                               << "\n";
1216           }
1217         }
1218       if (!corrupted)
1219         {
1220           // Apply the tail cancelation via the digital filter
1221           // (only for non-coorupted pads)
1222           DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1223         }
1224
1225       for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime))
1226         {
1227           // Store the amplitude of the digit if above threshold
1228           if (outADC[iTime] > 0)
1229             fDigits->SetData(iRow,iCol,iTime,TMath::Nint(outADC[iTime]));
1230           else
1231             fDigits->SetData(iRow,iCol,iTime,0);
1232         } // while itime
1233
1234     } // while irow icol
1235
1236   delete [] inADC;
1237   delete [] outADC;
1238
1239   return;
1240
1241 }
1242
1243 //_____________________________________________________________________________
1244 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1245                                   ,const Int_t n, const Int_t nexp)
1246 {
1247   //
1248   // Tail cancellation by deconvolution for PASA v4 TRF
1249   //
1250
1251   Double_t rates[2];
1252   Double_t coefficients[2];
1253
1254   // Initialization (coefficient = alpha, rates = lambda)
1255   Double_t r1 = 1.0;
1256   Double_t r2 = 1.0;
1257   Double_t c1 = 0.5;
1258   Double_t c2 = 0.5;
1259
1260   if (nexp == 1) {   // 1 Exponentials
1261     r1 = 1.156;
1262     r2 = 0.130;
1263     c1 = 0.066;
1264     c2 = 0.000;
1265   }
1266   if (nexp == 2) {   // 2 Exponentials
1267     Double_t par[4];
1268     fReconstructor->GetRecoParam()->GetTCParams(par);
1269     r1 = par[0];//1.156;
1270     r2 = par[1];//0.130;
1271     c1 = par[2];//0.114;
1272     c2 = par[3];//0.624;
1273   }
1274
1275   coefficients[0] = c1;
1276   coefficients[1] = c2;
1277
1278   Double_t dt = 0.1;
1279
1280   rates[0] = TMath::Exp(-dt/(r1));
1281   rates[1] = TMath::Exp(-dt/(r2));
1282
1283   Int_t i = 0;
1284   Int_t k = 0;
1285
1286   Double_t reminder[2];
1287   Double_t correction = 0.0;
1288   Double_t result     = 0.0;
1289
1290   // Attention: computation order is important
1291   for (k = 0; k < nexp; k++) {
1292     reminder[k] = 0.0;
1293   }
1294
1295   for (i = 0; i < n; i++) {
1296
1297     result    = (source[i] - correction);    // No rescaling
1298     target[i] = result;
1299
1300     for (k = 0; k < nexp; k++) {
1301       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1302     }
1303
1304     correction = 0.0;
1305     for (k = 0; k < nexp; k++) {
1306       correction += reminder[k];
1307     }
1308
1309   }
1310
1311 }
1312
1313 //_____________________________________________________________________________
1314 void AliTRDclusterizer::ResetRecPoints()
1315 {
1316   //
1317   // Resets the list of rec points
1318   //
1319
1320   if (fRecPoints) {
1321     fRecPoints->Delete();
1322     delete fRecPoints;
1323   }
1324 }
1325
1326 //_____________________________________________________________________________
1327 TClonesArray *AliTRDclusterizer::RecPoints()
1328 {
1329   //
1330   // Returns the list of rec points
1331   //
1332
1333   if (!fRecPoints) {
1334     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1335       // determine number of clusters which has to be allocated
1336       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1337
1338       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1339     }
1340     //SetClustersOwner(kTRUE);
1341     AliTRDReconstructor::SetClusters(0x0);
1342   }
1343   return fRecPoints;
1344
1345 }
1346
1347 //_____________________________________________________________________________
1348 TClonesArray *AliTRDclusterizer::TrackletsArray()
1349 {
1350   //
1351   // Returns the list of rec points
1352   //
1353
1354   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1355     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1356     //SetClustersOwner(kTRUE);
1357     //AliTRDReconstructor::SetTracklets(0x0);
1358   }
1359   return fTracklets;
1360
1361 }
1362