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