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