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