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