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