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