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