TRD/AliTRDclusterizer.cxx

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