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   if(fReconstructor->IsHLT()){
 670     fRawStream->SetSharedPadReadout(kFALSE);
 671     fRawStream->SetNoErrorWarning();
 672   }
 673
 674   AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
 675
 676   Int_t det    = 0;
 677   while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
 678     Bool_t iclusterBranch = kFALSE;
 679     if (fDigitsManager->GetIndexes(det)->HasEntry()){
 680       iclusterBranch = MakeClusters(det);
 681     }
 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   }
 702
 703   AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
 704   return kTRUE;
 705
 706 }
 707
 708 //_____________________________________________________________________________
 709 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 710 {
 711   //
 712   // Check if a pad is masked
 713   //
 714
 715   UChar_t status = 0;
 716
 717   if(signal>0 && TESTBIT(signal, 10)){
 718     CLRBIT(signal, 10);
 719     for(int ibit=0; ibit<4; ibit++){
 720       if(TESTBIT(signal, 11+ibit)){
 721         SETBIT(status, ibit);
 722         CLRBIT(signal, 11+ibit);
 723       }
 724     }
 725   }
 726   return status;
 727 }
 728
 729 //_____________________________________________________________________________
 730 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
 731   //
 732   // Set the pad status into out
 733   // First three bits are needed for the position encoding
 734   //
 735   out |= status << 3;
 736 }
 737
 738 //_____________________________________________________________________________
 739 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
 740   //
 741   // return the staus encoding of the corrupted pad
 742   //
 743   return static_cast<UChar_t>(encoding >> 3);
 744 }
 745
 746 //_____________________________________________________________________________
 747 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
 748   //
 749   // Return the position of the corruption
 750   //
 751   return encoding & 7;
 752 }
 753
 754 //_____________________________________________________________________________
 755 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 756 {
 757   //
 758   // Generates the cluster.
 759   //
 760
 761   // Get the digits
 762   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
 763   fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline();
 764
 765   // This is to take care of switched off super modules
 766   if (!fDigits->HasData()) return kFALSE;
 767
 768   fIndexes = fDigitsManager->GetIndexes(det);
 769   if (fIndexes->IsAllocated() == kFALSE) {
 770     AliError("Indexes do not exist!");
 771     return kFALSE;
 772   }
 773
 774   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
 775   if (!calibration) {
 776     AliFatal("No AliTRDcalibDB instance available\n");
 777     return kFALSE;
 778   }
 779
 780   if (!fReconstructor){
 781     AliError("Reconstructor not set\n");
 782     return kFALSE;
 783   }
 784
 785   fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
 786   fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
 787   fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
 788   fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
 789
 790   Int_t istack  = fIndexes->GetStack();
 791   fLayer  = fIndexes->GetLayer();
 792   Int_t isector = fIndexes->GetSM();
 793
 794   // Start clustering in the chamber
 795
 796   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 797   if (fDet != det) {
 798     AliError("Strange Detector number mismatch!");
 799     return kFALSE;
 800   }
 801
 802   AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
 803
 804   // TRD space point transformation
 805   fTransform->SetDetector(det);
 806
 807   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 808   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 809   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 810
 811   if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
 812     AddTrackletsToArray();
 813
 814   fColMax    = fDigits->GetNcol();
 815   //Int_t nRowMax    = fDigits->GetNrow();
 816   fTimeTotal = fDigits->GetNtime();
 817
 818   // Check consistency between OCDB and raw data
 819   if (fTimeTotal != calibration->GetNumberOfTimeBinsDCS()) {
 820     AliError(Form("Number of timebins does not match OCDB value (raw:%d, OCDB:%d)"
 821                  ,fTimeTotal,calibration->GetNumberOfTimeBinsDCS()));
 822   }
 823
 824   // Detector wise calibration object for the gain factors
 825   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 826   // Calibration object with pad wise values for the gain factors
 827   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 828   // Calibration value for chamber wise gain factor
 829   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 830
 831   // Detector wise calibration object for the noise
 832   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 833   // Calibration object with pad wise values for the noise
 834   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 835   // Calibration value for chamber wise noise
 836   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 837
 838   // Calibration object with the pad status
 839   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
 840
 841   SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
 842   SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
 843   SetBit(kHLT, fReconstructor->IsHLT());
 844
 845   firstClusterROC = -1;
 846   fClusterROC     =  0;
 847
 848   // Apply the gain and the tail cancelation via digital filter
 849   if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
 850
 851   MaxStruct curr, last;
 852   Int_t nMaximas = 0, nCorrupted = 0;
 853
 854   // Here the clusterfining is happening
 855
 856   for(curr.time = 0; curr.time < fTimeTotal; curr.time++){
 857     while(fIndexes->NextRCIndex(curr.row, curr.col)){
 858       if(IsMaximum(curr, curr.padStatus, &curr.signals[0])){
 859         if(last.row>-1){
 860           if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
 861           CreateCluster(last);
 862         }
 863         last=curr; curr.fivePad=kFALSE;
 864       }
 865     }
 866   }
 867   if(last.row>-1) CreateCluster(last);
 868
 869   if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
 870     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
 871     (*fDebugStream) << "MakeClusters"
 872       << "Detector="   << det
 873       << "NMaxima="    << nMaximas
 874       << "NClusters="  << fClusterROC
 875       << "NCorrupted=" << nCorrupted
 876       << "\n";
 877   }
 878   if (TestBit(kLabels)) AddLabels();
 879
 880   return kTRUE;
 881
 882 }
 883
 884 //_____________________________________________________________________________
 885 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
 886 {
 887   //
 888   // Returns true if this row,col,time combination is a maximum.
 889   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 890   //
 891
 892   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
 893   Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
 894   if(Signals[1] < fMaxThresh) return kFALSE;
 895
 896   Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
 897   if (Signals[1] < noiseMiddleThresh) return kFALSE;
 898
 899   if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
 900
 901   UChar_t status[3]={
 902     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
 903    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
 904    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
 905   };
 906
 907   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
 908   Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
 909   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
 910   Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
 911
 912   if(!(status[0] | status[1] | status[2])) {//all pads are good
 913     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 914       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
 915         if(Signals[0]<0)Signals[0]=0;
 916         if(Signals[2]<0)Signals[2]=0;
 917         Float_t  noiseSumThresh = fMinLeftRightCutSigma
 918           * fCalNoiseDetValue
 919           * fCalNoiseROC->GetValue(Max.col, Max.row);
 920         if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
 921         padStatus = 0;
 922         return kTRUE;
 923       }
 924     }
 925   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 926     if(Signals[0]<0)Signals[0]=0;
 927     if(Signals[2]<0)Signals[2]=0;
 928     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
 929       Signals[2]=0;
 930       SetPadStatus(status[2], padStatus);
 931       return kTRUE;
 932     }
 933     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
 934       Signals[0]=0;
 935       SetPadStatus(status[0], padStatus);
 936       return kTRUE;
 937     }
 938     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
 939       Signals[1] = (Short_t)(fMaxThresh + 0.5f);
 940       SetPadStatus(status[1], padStatus);
 941       return kTRUE;
 942     }
 943   }
 944   return kFALSE;
 945 }
 946
 947 //_____________________________________________________________________________
 948 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
 949 {
 950   //
 951   // Look for 5 pad cluster with minimum in the middle
 952   // Gives back the ratio
 953   //
 954
 955   if (ThisMax.col >= fColMax - 3) return kFALSE;
 956   Float_t gain;
 957   if (ThisMax.col < fColMax - 5){
 958     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
 959     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
 960       return kFALSE;
 961   }
 962   if (ThisMax.col > 1) {
 963     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
 964     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
 965       return kFALSE;
 966   }
 967
 968   const Float_t kEpsilon = 0.01;
 969   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
 970       NeighbourMax.signals[1], NeighbourMax.signals[2]};
 971
 972   // Unfold the two maxima and set the signal on
 973   // the overlapping pad to the ratio
 974   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
 975   ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
 976   NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
 977   ThisMax.fivePad=kTRUE;
 978   NeighbourMax.fivePad=kTRUE;
 979   return kTRUE;
 980
 981 }
 982
 983 //_____________________________________________________________________________
 984 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 985 {
 986   //
 987   // Creates a cluster at the given position and saves it in fRecPoints
 988   //
 989
 990   Int_t nPadCount = 1;
 991   Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
 992   if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
 993
 994   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
 995   cluster.SetNPads(nPadCount);
 996   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
 997   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
 998   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
 999
1000   cluster.SetFivePad(Max.fivePad);
1001   // set pads status for the cluster
1002   UChar_t maskPosition = GetCorruption(Max.padStatus);
1003   if (maskPosition) {
1004     cluster.SetPadMaskedPosition(maskPosition);
1005     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1006   }
1007
1008   // Transform the local cluster coordinates into calibrated
1009   // space point positions defined in the local tracking system.
1010   // Here the calibration for T0, Vdrift and ExB is applied as well.
1011   if(!fTransform->Transform(&cluster)) return;
1012   // Temporarily store the Max.Row, column and time bin of the center pad
1013   // Used to later on assign the track indices
1014   cluster.SetLabel(Max.row, 0);
1015   cluster.SetLabel(Max.col, 1);
1016   cluster.SetLabel(Max.time,2);
1017
1018   //needed for HLT reconstruction
1019   AddClusterToArray(&cluster);
1020
1021   // Store the index of the first cluster in the current ROC
1022   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1023
1024   fNoOfClusters++;
1025   fClusterROC++;
1026 }
1027
1028 //_____________________________________________________________________________
1029 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1030 {
1031   // Look to the right
1032   Int_t ii = 1;
1033   while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
1034     nPadCount++;
1035     ii++;
1036     if (Max.col < ii) break;
1037   }
1038   // Look to the left
1039   ii = 1;
1040   while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
1041     nPadCount++;
1042     ii++;
1043     if (Max.col+ii >= fColMax) break;
1044   }
1045
1046   // Store the amplitudes of the pads in the cluster for later analysis
1047   // and check whether one of these pads is masked in the database
1048   signals[2]=Max.signals[0];
1049   signals[3]=Max.signals[1];
1050   signals[4]=Max.signals[2];
1051   Float_t gain;
1052   for(Int_t i = 0; i<2; i++)
1053     {
1054       if(Max.col+i >= 3){
1055         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1056         signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1057       }
1058       if(Max.col+3-i < fColMax){
1059         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1060         signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1061       }
1062     }
1063   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1064     if ((jPad >= 0) && (jPad < fColMax))
1065       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1066       }*/
1067 }
1068
1069 //_____________________________________________________________________________
1070 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1071 {
1072   //
1073   // Add a cluster to the array
1074   //
1075
1076   Int_t n = RecPoints()->GetEntriesFast();
1077   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1078   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1079 }
1080
1081 //_____________________________________________________________________________
1082 void AliTRDclusterizer::AddTrackletsToArray()
1083 {
1084   //
1085   // Add the online tracklets of this chamber to the array
1086   //
1087
1088   UInt_t* trackletword;
1089   for(Int_t side=0; side<2; side++)
1090     {
1091       Int_t trkl=0;
1092       trackletword=fTrackletContainer[side];
1093       while(trackletword[trkl]>0){
1094         Int_t n = TrackletsArray()->GetEntriesFast();
1095         AliTRDtrackletWord tmp(trackletword[trkl]);
1096         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1097         trkl++;
1098       }
1099     }
1100 }
1101
1102 //_____________________________________________________________________________
1103 Bool_t AliTRDclusterizer::AddLabels()
1104 {
1105   //
1106   // Add the track indices to the found clusters
1107   //
1108
1109   const Int_t   kNclus  = 3;
1110   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1111   const Int_t   kNtrack = kNdict * kNclus;
1112
1113   Int_t iClusterROC = 0;
1114
1115   Int_t row  = 0;
1116   Int_t col  = 0;
1117   Int_t time = 0;
1118   Int_t iPad = 0;
1119
1120   // Temporary array to collect the track indices
1121   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1122
1123   // Loop through the dictionary arrays one-by-one
1124   // to keep memory consumption low
1125   AliTRDarrayDictionary *tracksIn = 0;  //mod
1126   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1127
1128     // tracksIn should be expanded beforehand!
1129     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1130
1131     // Loop though the clusters found in this ROC
1132     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1133
1134       AliTRDcluster *cluster = (AliTRDcluster *)
1135                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1136       row  = cluster->GetLabel(0);
1137       col  = cluster->GetLabel(1);
1138       time = cluster->GetLabel(2);
1139
1140       for (iPad = 0; iPad < kNclus; iPad++) {
1141         Int_t iPadCol = col - 1 + iPad;
1142         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1143         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1144       }
1145
1146     }
1147
1148   }
1149
1150   // Copy the track indices into the cluster
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     cluster->SetLabel(-9999,0);
1157     cluster->SetLabel(-9999,1);
1158     cluster->SetLabel(-9999,2);
1159
1160     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1161
1162   }
1163
1164   delete [] idxTracks;
1165
1166   return kTRUE;
1167
1168 }
1169
1170 //_____________________________________________________________________________
1171 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1172 {
1173   //
1174   // Method to unfold neighbouring maxima.
1175   // The charge ratio on the overlapping pad is calculated
1176   // until there is no more change within the range given by eps.
1177   // The resulting ratio is then returned to the calling method.
1178   //
1179
1180   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1181   if (!calibration) {
1182     AliError("No AliTRDcalibDB instance available\n");
1183     return kFALSE;
1184   }
1185
1186   Int_t   irc                = 0;
1187   Int_t   itStep             = 0;                 // Count iteration steps
1188
1189   Double_t ratio             = 0.5;               // Start value for ratio
1190   Double_t prevRatio         = 0.0;               // Store previous ratio
1191
1192   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1193   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1194   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1195
1196   // Start the iteration
1197   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1198
1199     itStep++;
1200     prevRatio = ratio;
1201
1202     // Cluster position according to charge ratio
1203     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1204                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1205     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1206                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1207
1208     // Set cluster charge ratio
1209     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1210     Double_t ampLeft  = padSignal[1] / newSignal[1];
1211     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1212     Double_t ampRight = padSignal[3] / newSignal[1];
1213
1214     // Apply pad response to parameters
1215     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1216     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1217
1218     // Calculate new overlapping ratio
1219     ratio = TMath::Min((Double_t) 1.0
1220                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1221
1222   }
1223
1224   return ratio;
1225
1226 }
1227
1228 //_____________________________________________________________________________
1229 void AliTRDclusterizer::TailCancelation()
1230 {
1231   //
1232   // Applies the tail cancelation and gain factors:
1233   // Transform fDigits to fDigits
1234   //
1235
1236   Int_t iRow  = 0;
1237   Int_t iCol  = 0;
1238   Int_t iTime = 0;
1239
1240   Double_t *inADC = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1241   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1242
1243   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1244   Bool_t debugStreaming = fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1245   while(fIndexes->NextRCIndex(iRow, iCol))
1246     {
1247       Bool_t corrupted = kFALSE;
1248       if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1249
1250       // Save data into the temporary processing array and substract the baseline,
1251       // since DeConvExp does not expect a baseline
1252       for (iTime = 0; iTime < fTimeTotal; iTime++)
1253         inADC[iTime]   = fDigits->GetData(iRow,iCol,iTime)-fBaseline;
1254
1255       if(debugStreaming){
1256         for (iTime = 0; iTime < fTimeTotal; iTime++)
1257           (*fDebugStream) << "TailCancellation"
1258                           << "col="  << iCol
1259                           << "row="  << iRow
1260                           << "time=" << iTime
1261                           << "inADC=" << inADC[iTime]
1262                           << "outADC=" << outADC[iTime]
1263                           << "corrupted=" << corrupted
1264                           << "\n";
1265       }
1266
1267       if (!corrupted)
1268         {
1269           // Apply the tail cancelation via the digital filter
1270           // (only for non-coorupted pads)
1271           DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1272         }
1273       else memcpy(&outADC[0],&inADC[0],fTimeTotal*sizeof(inADC[0]));
1274
1275       // Save tailcancalled data and add the baseline
1276       for(iTime = 0; iTime < fTimeTotal; iTime++)
1277         fDigits->SetData(iRow,iCol,iTime,(Short_t)(outADC[iTime] + fBaseline + 0.5));
1278
1279     } // while irow icol
1280
1281   delete [] inADC;
1282   delete [] outADC;
1283
1284   return;
1285
1286 }
1287
1288 //_____________________________________________________________________________
1289 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1290                                   ,const Int_t n, const Int_t nexp)
1291 {
1292   //
1293   // Tail cancellation by deconvolution for PASA v4 TRF
1294   //
1295
1296   Double_t rates[2];
1297   Double_t coefficients[2];
1298
1299   // Initialization (coefficient = alpha, rates = lambda)
1300   Double_t r1 = 1.0;
1301   Double_t r2 = 1.0;
1302   Double_t c1 = 0.5;
1303   Double_t c2 = 0.5;
1304
1305   if (nexp == 1) {   // 1 Exponentials
1306     r1 = 1.156;
1307     r2 = 0.130;
1308     c1 = 0.066;
1309     c2 = 0.000;
1310   }
1311   if (nexp == 2) {   // 2 Exponentials
1312     Double_t par[4];
1313     fReconstructor->GetRecoParam()->GetTCParams(par);
1314     r1 = par[0];//1.156;
1315     r2 = par[1];//0.130;
1316     c1 = par[2];//0.114;
1317     c2 = par[3];//0.624;
1318   }
1319
1320   coefficients[0] = c1;
1321   coefficients[1] = c2;
1322
1323   Double_t dt = 0.1;
1324
1325   rates[0] = TMath::Exp(-dt/(r1));
1326   rates[1] = TMath::Exp(-dt/(r2));
1327
1328   Int_t i = 0;
1329   Int_t k = 0;
1330
1331   Double_t reminder[2];
1332   Double_t correction = 0.0;
1333   Double_t result     = 0.0;
1334
1335   // Attention: computation order is important
1336   for (k = 0; k < nexp; k++) {
1337     reminder[k] = 0.0;
1338   }
1339
1340   for (i = 0; i < n; i++) {
1341
1342     result    = (source[i] - correction);    // No rescaling
1343     target[i] = result;
1344
1345     for (k = 0; k < nexp; k++) {
1346       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1347     }
1348
1349     correction = 0.0;
1350     for (k = 0; k < nexp; k++) {
1351       correction += reminder[k];
1352     }
1353
1354   }
1355
1356 }
1357
1358 //_____________________________________________________________________________
1359 void AliTRDclusterizer::ResetRecPoints()
1360 {
1361   //
1362   // Resets the list of rec points
1363   //
1364
1365   if (fRecPoints) {
1366     fRecPoints->Delete();
1367     delete fRecPoints;
1368   }
1369 }
1370
1371 //_____________________________________________________________________________
1372 TClonesArray *AliTRDclusterizer::RecPoints()
1373 {
1374   //
1375   // Returns the list of rec points
1376   //
1377
1378   if (!fRecPoints) {
1379     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1380       // determine number of clusters which has to be allocated
1381       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1382
1383       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1384     }
1385     //SetClustersOwner(kTRUE);
1386     AliTRDReconstructor::SetClusters(0x0);
1387   }
1388   return fRecPoints;
1389
1390 }
1391
1392 //_____________________________________________________________________________
1393 TClonesArray *AliTRDclusterizer::TrackletsArray()
1394 {
1395   //
1396   // Returns the list of rec points
1397   //
1398
1399   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1400     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1401     //SetClustersOwner(kTRUE);
1402     //AliTRDReconstructor::SetTracklets(0x0);
1403   }
1404   return fTracklets;
1405
1406 }
1407