TRD/AliTRDclusterizer.cxx

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