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