TRD/AliTRDclusterizer.cxx

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