TRD/AliTRDclusterizer.cxx

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