TRD/AliTRDclusterizer.cxx

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