TRD/AliTRDclusterizer.cxx

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