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