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     ((AliTRDrawStream*)fRawStream)->DisableErrorStorage();
 645   }
 646
 647   AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
 648
 649   UInt_t det = 0;
 650   while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) < AliTRDgeometry::kNdet){
 651     if (fDigitsManager->GetIndexes(det)->HasEntry()){
 652       MakeClusters(det);
 653       fDigitsManager->ClearArrays(det);
 654     }
 655
 656     if (!fReconstructor->IsWritingTracklets()) continue;
 657     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
 658   }
 659
 660   if (fReconstructor->IsWritingTracklets()) {
 661     if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
 662       if (trklLoader) {
 663         if (AliTreeLoader *trklTreeLoader = (AliTreeLoader*) trklLoader->GetBaseLoader("tracklets-raw"))
 664           trklTreeLoader->WriteData("OVERWRITE");
 665         trklLoader->UnloadAll();
 666       }
 667     }
 668   }
 669
 670   if (fTrackletContainer){
 671     delete [] fTrackletContainer[0];
 672     delete [] fTrackletContainer[1];
 673     delete [] fTrackletContainer;
 674     fTrackletContainer = NULL;
 675   }
 676
 677   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 678
 679   if(!TestBit(knewDM)){
 680     delete fDigitsManager;
 681     fDigitsManager = NULL;
 682     delete fRawStream;
 683     fRawStream = NULL;
 684   }
 685
 686   AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
 687   return kTRUE;
 688
 689 }
 690
 691 //_____________________________________________________________________________
 692 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 693 {
 694   //
 695   // Check if a pad is masked
 696   //
 697
 698   UChar_t status = 0;
 699
 700   if(signal>0 && TESTBIT(signal, 10)){
 701     CLRBIT(signal, 10);
 702     for(int ibit=0; ibit<4; ibit++){
 703       if(TESTBIT(signal, 11+ibit)){
 704         SETBIT(status, ibit);
 705         CLRBIT(signal, 11+ibit);
 706       }
 707     }
 708   }
 709   return status;
 710 }
 711
 712 //_____________________________________________________________________________
 713 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
 714   //
 715   // Set the pad status into out
 716   // First three bits are needed for the position encoding
 717   //
 718   out |= status << 3;
 719 }
 720
 721 //_____________________________________________________________________________
 722 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
 723   //
 724   // return the staus encoding of the corrupted pad
 725   //
 726   return static_cast<UChar_t>(encoding >> 3);
 727 }
 728
 729 //_____________________________________________________________________________
 730 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
 731   //
 732   // Return the position of the corruption
 733   //
 734   return encoding & 7;
 735 }
 736
 737 //_____________________________________________________________________________
 738 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 739 {
 740   //
 741   // Generates the cluster.
 742   //
 743
 744   // Get the digits
 745   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
 746   fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline(det);
 747
 748   // This is to take care of switched off super modules
 749   if (!fDigits->HasData()) return kFALSE;
 750
 751   fIndexes = fDigitsManager->GetIndexes(det);
 752   if (fIndexes->IsAllocated() == kFALSE) {
 753     AliError("Indexes do not exist!");
 754     return kFALSE;
 755   }
 756
 757   AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance();
 758   if (!calibration) {
 759     AliFatal("No AliTRDcalibDB instance available\n");
 760     return kFALSE;
 761   }
 762
 763   if (!fReconstructor){
 764     AliError("Reconstructor not set\n");
 765     return kFALSE;
 766   }
 767
 768   const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
 769
 770   fMaxThresh            = (Short_t)recoParam->GetClusMaxThresh();
 771   fMaxThreshTest        = (Short_t)(recoParam->GetClusMaxThresh()/2+fBaseline);
 772   fSigThresh            = (Short_t)recoParam->GetClusSigThresh();
 773   fMinMaxCutSigma       = recoParam->GetMinMaxCutSigma();
 774   fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma();
 775   const Int_t iEveryNTB = recoParam->GetRecEveryNTB();
 776
 777   Int_t istack  = fIndexes->GetStack();
 778   fLayer  = fIndexes->GetLayer();
 779   Int_t isector = fIndexes->GetSM();
 780
 781   // Start clustering in the chamber
 782
 783   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 784   if (fDet != det) {
 785     AliError("Strange Detector number Missmatch!");
 786     return kFALSE;
 787   }
 788
 789   AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
 790
 791   // TRD space point transformation
 792   fTransform->SetDetector(det);
 793
 794   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 795   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 796   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 797
 798   if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
 799     AddTrackletsToArray();
 800
 801   fColMax    = fDigits->GetNcol();
 802   fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
 803
 804   // Check consistency between OCDB and raw data
 805   Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
 806   if(fReconstructor->IsHLT()){
 807     if((nTimeOCDB > -1) && (fTimeTotal != nTimeOCDB)){
 808       AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value"
 809                       ,fTimeTotal,nTimeOCDB));
 810     }
 811   }else{
 812     if(nTimeOCDB == -1){
 813       AliWarning("Undefined number of timebins in OCDB, using value from raw data.");
 814       if(!fTimeTotal>0){
 815         AliError("Number of timebins in raw data is negative, skipping chamber!");
 816         return kFALSE;
 817       }
 818     }else if(nTimeOCDB == -2){
 819       AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!");
 820       return kFALSE;
 821     }else if(fTimeTotal != nTimeOCDB){
 822       AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber!"
 823                     ,fTimeTotal,nTimeOCDB));
 824       return kFALSE;
 825     }
 826   }
 827
 828   // Detector wise calibration object for the gain factors
 829   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 830   // Calibration object with pad wise values for the gain factors
 831   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 832   // Calibration value for chamber wise gain factor
 833   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 834
 835   // Detector wise calibration object for the noise
 836   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 837   // Calibration object with pad wise values for the noise
 838   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 839   // Calibration value for chamber wise noise
 840   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 841
 842   // Calibration object with the pad status
 843   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
 844
 845   firstClusterROC = -1;
 846   fClusterROC     =  0;
 847
 848   SetBit(kLUT, recoParam->UseLUT());
 849   SetBit(kGAUS, recoParam->UseGAUS());
 850
 851   // Apply the gain and the tail cancelation via digital filter
 852   // Use the configuration from the DCS to find out whether online
 853   // tail cancellation was applied
 854   if(!calibration->HasOnlineTailCancellation()) TailCancelation(recoParam);
 855
 856   MaxStruct curr, last;
 857   Int_t nMaximas = 0, nCorrupted = 0;
 858
 859   // Here the clusterfining is happening
 860
 861   for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
 862     while(fIndexes->NextRCIndex(curr.row, curr.col)){
 863       if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
 864         if(last.row>-1){
 865           if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
 866           CreateCluster(last);
 867         }
 868         last=curr; curr.fivePad=kFALSE;
 869       }
 870     }
 871   }
 872   if(last.row>-1) CreateCluster(last);
 873
 874   if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
 875     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
 876     (*fDebugStream) << "MakeClusters"
 877       << "Detector="   << det
 878       << "NMaxima="    << nMaximas
 879       << "NClusters="  << fClusterROC
 880       << "NCorrupted=" << nCorrupted
 881       << "\n";
 882   }
 883   if (TestBit(kLabels)) AddLabels();
 884
 885   return kTRUE;
 886
 887 }
 888
 889 //_____________________________________________________________________________
 890 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
 891 {
 892   //
 893   // Returns true if this row,col,time combination is a maximum.
 894   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 895   //
 896
 897   Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
 898   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
 899   tmp = (fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f;
 900   Signals[1] = (Short_t)TMath::Min(tmp, kMaxShortVal);
 901   if(Signals[1] <= fMaxThresh) return kFALSE;
 902
 903   if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
 904
 905   Short_t noiseMiddleThresh = (Short_t)(fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row));
 906   if (Signals[1] <= noiseMiddleThresh) return kFALSE;
 907
 908   UChar_t status[3]={
 909     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
 910    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
 911    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
 912   };
 913
 914   Short_t signal(0);
 915   if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
 916     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
 917     tmp = (signal - fBaseline) / gain + 0.5f;
 918     Signals[0] = (Short_t)TMath::Min(tmp, kMaxShortVal);
 919   } else Signals[0] = 0;
 920   if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
 921     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
 922     tmp = (signal - fBaseline) / gain + 0.5f;
 923     Signals[2] = (Short_t)TMath::Min(tmp, kMaxShortVal);
 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 tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
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     tmp = (signal - fBaseline) / gain + 0.5f;
1060     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1061     if(signal<fSigThresh) break; // signal under threshold
1062     nPadCount++;
1063     if(ipad<=3) signals[3 - ipad] = signal;
1064     ipad++;
1065   }
1066   ipad=1;
1067   // Look to the left
1068   while((jpad = Max.col+ipad)<fColMax){
1069     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1070     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1071     tmp = (signal - fBaseline) / gain + 0.5f;
1072     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1073     if(signal<fSigThresh) break; // signal under threshold
1074     nPadCount++;
1075     if(ipad<=3) signals[3 + ipad] = signal;
1076     ipad++;
1077   }
1078
1079   AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1080     , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
1081 }
1082
1083 //_____________________________________________________________________________
1084 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1085 {
1086   //
1087   // Add a cluster to the array
1088   //
1089
1090   Int_t n = RecPoints()->GetEntriesFast();
1091   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1092   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1093 }
1094
1095 //_____________________________________________________________________________
1096 void AliTRDclusterizer::AddTrackletsToArray()
1097 {
1098   //
1099   // Add the online tracklets of this chamber to the array
1100   //
1101
1102   UInt_t* trackletword;
1103   for(Int_t side=0; side<2; side++)
1104     {
1105       Int_t trkl=0;
1106       trackletword=fTrackletContainer[side];
1107       while(trackletword[trkl]>0){
1108         Int_t n = TrackletsArray()->GetEntriesFast();
1109         AliTRDtrackletWord tmp(trackletword[trkl]);
1110         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1111         trkl++;
1112       }
1113     }
1114 }
1115
1116 //_____________________________________________________________________________
1117 Bool_t AliTRDclusterizer::AddLabels()
1118 {
1119   //
1120   // Add the track indices to the found clusters
1121   //
1122
1123   const Int_t   kNclus  = 3;
1124   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1125   const Int_t   kNtrack = kNdict * kNclus;
1126
1127   Int_t iClusterROC = 0;
1128
1129   Int_t row  = 0;
1130   Int_t col  = 0;
1131   Int_t time = 0;
1132   Int_t iPad = 0;
1133
1134   // Temporary array to collect the track indices
1135   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1136
1137   // Loop through the dictionary arrays one-by-one
1138   // to keep memory consumption low
1139   AliTRDarrayDictionary *tracksIn = 0;  //mod
1140   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1141
1142     // tracksIn should be expanded beforehand!
1143     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1144
1145     // Loop though the clusters found in this ROC
1146     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1147
1148       AliTRDcluster *cluster = (AliTRDcluster *)
1149                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1150       row  = cluster->GetLabel(0);
1151       col  = cluster->GetLabel(1);
1152       time = cluster->GetLabel(2);
1153
1154       for (iPad = 0; iPad < kNclus; iPad++) {
1155         Int_t iPadCol = col - 1 + iPad;
1156         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1157         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1158       }
1159
1160     }
1161
1162   }
1163
1164   // Copy the track indices into the cluster
1165   // Loop though the clusters found in this ROC
1166   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1167
1168     AliTRDcluster *cluster = (AliTRDcluster *)
1169       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1170     cluster->SetLabel(-9999,0);
1171     cluster->SetLabel(-9999,1);
1172     cluster->SetLabel(-9999,2);
1173
1174     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1175
1176   }
1177
1178   delete [] idxTracks;
1179
1180   return kTRUE;
1181
1182 }
1183
1184 //_____________________________________________________________________________
1185 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1186 {
1187   //
1188   // Method to unfold neighbouring maxima.
1189   // The charge ratio on the overlapping pad is calculated
1190   // until there is no more change within the range given by eps.
1191   // The resulting ratio is then returned to the calling method.
1192   //
1193
1194   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1195   if (!calibration) {
1196     AliError("No AliTRDcalibDB instance available\n");
1197     return kFALSE;
1198   }
1199
1200   Int_t   irc                = 0;
1201   Int_t   itStep             = 0;                 // Count iteration steps
1202
1203   Double_t ratio             = 0.5;               // Start value for ratio
1204   Double_t prevRatio         = 0.0;               // Store previous ratio
1205
1206   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1207   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1208   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1209
1210   // Start the iteration
1211   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1212
1213     itStep++;
1214     prevRatio = ratio;
1215
1216     // Cluster position according to charge ratio
1217     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1218                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1219     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1220                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1221
1222     // Set cluster charge ratio
1223     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1224     Double_t ampLeft  = padSignal[1] / newSignal[1];
1225     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1226     Double_t ampRight = padSignal[3] / newSignal[1];
1227
1228     // Apply pad response to parameters
1229     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1230     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1231
1232     // Calculate new overlapping ratio
1233     ratio = TMath::Min((Double_t) 1.0
1234                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1235
1236   }
1237
1238   return ratio;
1239
1240 }
1241
1242 //_____________________________________________________________________________
1243 void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1244 {
1245   //
1246   // Applies the tail cancelation
1247   //
1248
1249   Int_t nexp = recoParam->GetTCnexp();
1250   if(!nexp) return;
1251
1252   Int_t iRow  = 0;
1253   Int_t iCol  = 0;
1254   Int_t iTime = 0;
1255
1256   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1257   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1258   while(fIndexes->NextRCIndex(iRow, iCol))
1259     {
1260       // if corrupted then don't make the tail cancallation
1261       if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1262
1263       if(debugStreaming){
1264         for (iTime = 0; iTime < fTimeTotal; iTime++)
1265           (*fDebugStream) << "TailCancellation"
1266                           << "col="  << iCol
1267                           << "row="  << iRow
1268                           << "\n";
1269       }
1270
1271       // Apply the tail cancelation via the digital filter
1272       DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1273
1274     } // while irow icol
1275
1276   return;
1277
1278 }
1279
1280 //_____________________________________________________________________________
1281 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1282 {
1283   //
1284   // Tail cancellation by deconvolution for PASA v4 TRF
1285   //
1286
1287   Float_t rates[2];
1288   Float_t coefficients[2];
1289
1290   // Initialization (coefficient = alpha, rates = lambda)
1291   Float_t r1 = 1.0;
1292   Float_t r2 = 1.0;
1293   Float_t c1 = 0.5;
1294   Float_t c2 = 0.5;
1295
1296   if (nexp == 1) {   // 1 Exponentials
1297     r1 = 1.156;
1298     r2 = 0.130;
1299     c1 = 0.066;
1300     c2 = 0.000;
1301   }
1302   if (nexp == 2) {   // 2 Exponentials
1303     Double_t par[4];
1304     fReconstructor->GetRecoParam()->GetTCParams(par);
1305     r1 = par[0];//1.156;
1306     r2 = par[1];//0.130;
1307     c1 = par[2];//0.114;
1308     c2 = par[3];//0.624;
1309   }
1310
1311   coefficients[0] = c1;
1312   coefficients[1] = c2;
1313
1314   Double_t dt = 0.1;
1315
1316   rates[0] = TMath::Exp(-dt/(r1));
1317   rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
1318
1319   Float_t reminder[2] = { .0, .0 };
1320   Float_t correction = 0.0;
1321   Float_t result     = 0.0;
1322
1323   for (int i = 0; i < nTime; i++) {
1324
1325     result = arr[i] - correction - fBaseline;    // No rescaling
1326     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1327
1328     correction = 0.0;
1329     for (int k = 0; k < 2; k++) {
1330       correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1331     }
1332
1333   }
1334
1335 }
1336
1337 //_____________________________________________________________________________
1338 void AliTRDclusterizer::ResetRecPoints()
1339 {
1340   //
1341   // Resets the list of rec points
1342   //
1343
1344   if (fRecPoints) {
1345     fRecPoints->Clear();
1346     fNoOfClusters = 0;
1347     //    delete fRecPoints;
1348   }
1349 }
1350
1351 //_____________________________________________________________________________
1352 TClonesArray *AliTRDclusterizer::RecPoints()
1353 {
1354   //
1355   // Returns the list of rec points
1356   //
1357
1358   if (!fRecPoints) {
1359     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1360       // determine number of clusters which has to be allocated
1361       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1362
1363       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1364     }
1365     //SetClustersOwner(kTRUE);
1366     AliTRDReconstructor::SetClusters(0x0);
1367   }
1368   return fRecPoints;
1369
1370 }
1371
1372 //_____________________________________________________________________________
1373 TClonesArray *AliTRDclusterizer::TrackletsArray()
1374 {
1375   //
1376   // Returns the list of rec points
1377   //
1378
1379   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1380     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1381     //SetClustersOwner(kTRUE);
1382     //AliTRDReconstructor::SetTracklets(0x0);
1383   }
1384   return fTracklets;
1385
1386 }
1387