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