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