TRD/AliTRDclusterizer.cxx

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