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     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   Short_t signal(0);
 917   if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
 918     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
 919     Signals[0] = (Short_t)((signal - fBaseline) / gain + 0.5f);
 920   } else Signals[0] = 0;
 921   if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
 922     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
 923     Signals[2] = (Short_t)((signal - fBaseline) / gain + 0.5f);
 924   } else Signals[2] = 0;
 925
 926   if(!(status[0] | status[1] | status[2])) {//all pads are good
 927     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 928       if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
 929         if(Signals[0]<0)Signals[0]=0;
 930         if(Signals[2]<0)Signals[2]=0;
 931         Short_t noiseSumThresh = (Short_t)(fMinLeftRightCutSigma * fCalNoiseDetValue
 932                                            * fCalNoiseROC->GetValue(Max.col, Max.row));
 933         if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
 934         padStatus = 0;
 935         return kTRUE;
 936       }
 937     }
 938   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 939     if(Signals[0]<0)Signals[0]=0;
 940     if(Signals[2]<0)Signals[2]=0;
 941     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
 942       Signals[2]=0;
 943       SetPadStatus(status[2], padStatus);
 944       return kTRUE;
 945     }
 946     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
 947       Signals[0]=0;
 948       SetPadStatus(status[0], padStatus);
 949       return kTRUE;
 950     }
 951     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
 952       Signals[1] = fMaxThresh;
 953       SetPadStatus(status[1], padStatus);
 954       return kTRUE;
 955     }
 956   }
 957   return kFALSE;
 958 }
 959
 960 //_____________________________________________________________________________
 961 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
 962 {
 963   //
 964   // Look for 5 pad cluster with minimum in the middle
 965   // Gives back the ratio
 966   //
 967
 968   if (ThisMax.col >= fColMax - 3) return kFALSE;
 969   Float_t gain;
 970   if (ThisMax.col < fColMax - 5){
 971     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
 972     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
 973       return kFALSE;
 974   }
 975   if (ThisMax.col > 1) {
 976     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
 977     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
 978       return kFALSE;
 979   }
 980
 981   const Float_t kEpsilon = 0.01;
 982   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
 983       NeighbourMax.signals[1], NeighbourMax.signals[2]};
 984
 985   // Unfold the two maxima and set the signal on
 986   // the overlapping pad to the ratio
 987   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
 988   ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
 989   NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
 990   ThisMax.fivePad=kTRUE;
 991   NeighbourMax.fivePad=kTRUE;
 992   return kTRUE;
 993
 994 }
 995
 996 //_____________________________________________________________________________
 997 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 998 {
 999   //
1000   // Creates a cluster at the given position and saves it in fRecPoints
1001   //
1002
1003   Int_t nPadCount = 1;
1004   Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
1005   if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
1006
1007   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
1008   cluster.SetNPads(nPadCount);
1009   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1010   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1011   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
1012
1013   cluster.SetFivePad(Max.fivePad);
1014   // set pads status for the cluster
1015   UChar_t maskPosition = GetCorruption(Max.padStatus);
1016   if (maskPosition) {
1017     cluster.SetPadMaskedPosition(maskPosition);
1018     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1019   }
1020   cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
1021
1022   // Transform the local cluster coordinates into calibrated
1023   // space point positions defined in the local tracking system.
1024   // Here the calibration for T0, Vdrift and ExB is applied as well.
1025   if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
1026
1027   // Temporarily store the Max.Row, column and time bin of the center pad
1028   // Used to later on assign the track indices
1029   cluster.SetLabel(Max.row, 0);
1030   cluster.SetLabel(Max.col, 1);
1031   cluster.SetLabel(Max.time,2);
1032
1033   //needed for HLT reconstruction
1034   AddClusterToArray(&cluster);
1035
1036   // Store the index of the first cluster in the current ROC
1037   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1038
1039   fNoOfClusters++;
1040   fClusterROC++;
1041 }
1042
1043 //_____________________________________________________________________________
1044 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1045 {
1046 // Calculate number of pads/cluster and
1047 // ADC signals at position 0, 1, 5 and 6
1048
1049   Float_t gain(1.); Short_t signal(0.);
1050   // Store the amplitudes of the pads in the cluster for later analysis
1051   // and check whether one of these pads is masked in the database
1052   signals[3]=Max.signals[1];
1053   Int_t ipad(1), jpad(0);
1054   // Look to the right
1055   while((jpad = Max.col-ipad)){
1056     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1057     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1058     signal = (Short_t)((signal - fBaseline) / gain + 0.5f);
1059     if(signal<fSigThresh) break; // signal under threshold
1060     nPadCount++;
1061     if(ipad<=3) signals[3 - ipad] = signal;
1062     ipad++;
1063   }
1064   ipad=1;
1065   // Look to the left
1066   while((jpad = Max.col+ipad)<fColMax){
1067     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1068     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1069     signal = (Short_t)((signal - fBaseline) / gain + 0.5f);
1070     if(signal<fSigThresh) break; // signal under threshold
1071     nPadCount++;
1072     if(ipad<=3) signals[3 + ipad] = signal;
1073     ipad++;
1074   }
1075
1076   AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1077     , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
1078 }
1079
1080 //_____________________________________________________________________________
1081 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1082 {
1083   //
1084   // Add a cluster to the array
1085   //
1086
1087   Int_t n = RecPoints()->GetEntriesFast();
1088   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1089   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1090 }
1091
1092 //_____________________________________________________________________________
1093 void AliTRDclusterizer::AddTrackletsToArray()
1094 {
1095   //
1096   // Add the online tracklets of this chamber to the array
1097   //
1098
1099   UInt_t* trackletword;
1100   for(Int_t side=0; side<2; side++)
1101     {
1102       Int_t trkl=0;
1103       trackletword=fTrackletContainer[side];
1104       while(trackletword[trkl]>0){
1105         Int_t n = TrackletsArray()->GetEntriesFast();
1106         AliTRDtrackletWord tmp(trackletword[trkl]);
1107         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1108         trkl++;
1109       }
1110     }
1111 }
1112
1113 //_____________________________________________________________________________
1114 Bool_t AliTRDclusterizer::AddLabels()
1115 {
1116   //
1117   // Add the track indices to the found clusters
1118   //
1119
1120   const Int_t   kNclus  = 3;
1121   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1122   const Int_t   kNtrack = kNdict * kNclus;
1123
1124   Int_t iClusterROC = 0;
1125
1126   Int_t row  = 0;
1127   Int_t col  = 0;
1128   Int_t time = 0;
1129   Int_t iPad = 0;
1130
1131   // Temporary array to collect the track indices
1132   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1133
1134   // Loop through the dictionary arrays one-by-one
1135   // to keep memory consumption low
1136   AliTRDarrayDictionary *tracksIn = 0;  //mod
1137   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1138
1139     // tracksIn should be expanded beforehand!
1140     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1141
1142     // Loop though the clusters found in this ROC
1143     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1144
1145       AliTRDcluster *cluster = (AliTRDcluster *)
1146                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1147       row  = cluster->GetLabel(0);
1148       col  = cluster->GetLabel(1);
1149       time = cluster->GetLabel(2);
1150
1151       for (iPad = 0; iPad < kNclus; iPad++) {
1152         Int_t iPadCol = col - 1 + iPad;
1153         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1154         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1155       }
1156
1157     }
1158
1159   }
1160
1161   // Copy the track indices into the cluster
1162   // Loop though the clusters found in this ROC
1163   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1164
1165     AliTRDcluster *cluster = (AliTRDcluster *)
1166       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1167     cluster->SetLabel(-9999,0);
1168     cluster->SetLabel(-9999,1);
1169     cluster->SetLabel(-9999,2);
1170
1171     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1172
1173   }
1174
1175   delete [] idxTracks;
1176
1177   return kTRUE;
1178
1179 }
1180
1181 //_____________________________________________________________________________
1182 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1183 {
1184   //
1185   // Method to unfold neighbouring maxima.
1186   // The charge ratio on the overlapping pad is calculated
1187   // until there is no more change within the range given by eps.
1188   // The resulting ratio is then returned to the calling method.
1189   //
1190
1191   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1192   if (!calibration) {
1193     AliError("No AliTRDcalibDB instance available\n");
1194     return kFALSE;
1195   }
1196
1197   Int_t   irc                = 0;
1198   Int_t   itStep             = 0;                 // Count iteration steps
1199
1200   Double_t ratio             = 0.5;               // Start value for ratio
1201   Double_t prevRatio         = 0.0;               // Store previous ratio
1202
1203   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1204   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1205   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1206
1207   // Start the iteration
1208   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1209
1210     itStep++;
1211     prevRatio = ratio;
1212
1213     // Cluster position according to charge ratio
1214     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1215                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1216     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1217                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1218
1219     // Set cluster charge ratio
1220     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1221     Double_t ampLeft  = padSignal[1] / newSignal[1];
1222     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1223     Double_t ampRight = padSignal[3] / newSignal[1];
1224
1225     // Apply pad response to parameters
1226     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1227     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1228
1229     // Calculate new overlapping ratio
1230     ratio = TMath::Min((Double_t) 1.0
1231                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1232
1233   }
1234
1235   return ratio;
1236
1237 }
1238
1239 //_____________________________________________________________________________
1240 void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1241 {
1242   //
1243   // Applies the tail cancelation
1244   //
1245
1246   Int_t iRow  = 0;
1247   Int_t iCol  = 0;
1248   Int_t iTime = 0;
1249
1250   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1251   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1252   Int_t nexp = recoParam->GetTCnexp();
1253   while(fIndexes->NextRCIndex(iRow, iCol))
1254     {
1255       // if corrupted then don't make the tail cancallation
1256       if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1257
1258       if(debugStreaming){
1259         for (iTime = 0; iTime < fTimeTotal; iTime++)
1260           (*fDebugStream) << "TailCancellation"
1261                           << "col="  << iCol
1262                           << "row="  << iRow
1263                           << "\n";
1264       }
1265
1266       // Apply the tail cancelation via the digital filter
1267       DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1268
1269     } // while irow icol
1270
1271   return;
1272
1273 }
1274
1275 //_____________________________________________________________________________
1276 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1277 {
1278   //
1279   // Tail cancellation by deconvolution for PASA v4 TRF
1280   //
1281
1282   Float_t rates[2];
1283   Float_t coefficients[2];
1284
1285   // Initialization (coefficient = alpha, rates = lambda)
1286   Float_t r1 = 1.0;
1287   Float_t r2 = 1.0;
1288   Float_t c1 = 0.5;
1289   Float_t c2 = 0.5;
1290
1291   if (nexp == 1) {   // 1 Exponentials
1292     r1 = 1.156;
1293     r2 = 0.130;
1294     c1 = 0.066;
1295     c2 = 0.000;
1296   }
1297   if (nexp == 2) {   // 2 Exponentials
1298     Double_t par[4];
1299     fReconstructor->GetRecoParam()->GetTCParams(par);
1300     r1 = par[0];//1.156;
1301     r2 = par[1];//0.130;
1302     c1 = par[2];//0.114;
1303     c2 = par[3];//0.624;
1304   }
1305
1306   coefficients[0] = c1;
1307   coefficients[1] = c2;
1308
1309   Double_t dt = 0.1;
1310
1311   rates[0] = TMath::Exp(-dt/(r1));
1312   rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
1313
1314   Float_t reminder[2] = { .0, .0 };
1315   Float_t correction = 0.0;
1316   Float_t result     = 0.0;
1317
1318   for (int i = 0; i < nTime; i++) {
1319
1320     result = arr[i] - correction - fBaseline;    // No rescaling
1321     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1322
1323     correction = 0.0;
1324     for (int k = 0; k < 2; k++) {
1325       correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1326     }
1327
1328   }
1329
1330 }
1331
1332 //_____________________________________________________________________________
1333 void AliTRDclusterizer::ResetRecPoints()
1334 {
1335   //
1336   // Resets the list of rec points
1337   //
1338
1339   if (fRecPoints) {
1340     fRecPoints->Clear();
1341     fNoOfClusters = 0;
1342     //    delete fRecPoints;
1343   }
1344 }
1345
1346 //_____________________________________________________________________________
1347 TClonesArray *AliTRDclusterizer::RecPoints()
1348 {
1349   //
1350   // Returns the list of rec points
1351   //
1352
1353   if (!fRecPoints) {
1354     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1355       // determine number of clusters which has to be allocated
1356       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1357
1358       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1359     }
1360     //SetClustersOwner(kTRUE);
1361     AliTRDReconstructor::SetClusters(0x0);
1362   }
1363   return fRecPoints;
1364
1365 }
1366
1367 //_____________________________________________________________________________
1368 TClonesArray *AliTRDclusterizer::TrackletsArray()
1369 {
1370   //
1371   // Returns the list of rec points
1372   //
1373
1374   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1375     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1376     //SetClustersOwner(kTRUE);
1377     //AliTRDReconstructor::SetTracklets(0x0);
1378   }
1379   return fTracklets;
1380
1381 }
1382