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