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