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
 547   AliDataLoader *trackLoader = loader->GetDataLoader("gtutracks");
 548   if (!trackLoader) {
 549       return kFALSE;
 550   }
 551
 552   trackLoader->Load();
 553
 554   TTree *trackTree = trackLoader->Tree();
 555   if (!trackTree) {
 556     return kFALSE;
 557   }
 558
 559   TClonesArray *trackArray = TracksArray();
 560   AliTRDtrackGTU *trk = 0x0;
 561   trackTree->SetBranchAddress("TRDtrackGTU", &trk);
 562   for (Int_t iTrack = 0; iTrack < trackTree->GetEntries(); iTrack++) {
 563     trackTree->GetEntry(iTrack);
 564     new ((*trackArray)[trackArray->GetEntries()]) AliESDTrdTrack(*(trk->CreateTrdTrack()));
 565   }
 566
 567   return kTRUE;
 568 }
 569
 570 //_____________________________________________________________________________
 571 Bool_t AliTRDclusterizer::MakeClusters()
 572 {
 573   //
 574   // Creates clusters from digits
 575   //
 576
 577   // Propagate info from the digits manager
 578   if (TestBit(kLabels)){
 579     SetBit(kLabels, fDigitsManager->UsesDictionaries());
 580   }
 581
 582   Bool_t fReturn = kTRUE;
 583   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
 584
 585     AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
 586     // This is to take care of switched off super modules
 587     if (!digitsIn->HasData()) continue;
 588     digitsIn->Expand();
 589     digitsIn->DeleteNegatives();  // Restore digits array to >=0 values
 590     AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 591     if (indexes->IsAllocated() == kFALSE){
 592       fDigitsManager->BuildIndexes(i);
 593     }
 594
 595     Bool_t fR = kFALSE;
 596     if (indexes->HasEntry()){
 597       if (TestBit(kLabels)){
 598         for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
 599           AliTRDarrayDictionary *tracksIn = 0; //mod
 600           tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict);  //mod
 601           tracksIn->Expand();
 602         }
 603       }
 604       fR = MakeClusters(i);
 605       fReturn = fR && fReturn;
 606     }
 607
 608     //if (fR == kFALSE){
 609     //  if(IsWritingClusters()) WriteClusters(i);
 610     //  ResetRecPoints();
 611     //}
 612
 613     // Clear arrays of this chamber, to prepare for next event
 614     fDigitsManager->ClearArrays(i);
 615   }
 616
 617   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 618
 619   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
 620
 621   return fReturn;
 622
 623 }
 624
 625 //_____________________________________________________________________________
 626 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
 627 {
 628   //
 629   // Creates clusters from raw data
 630   //
 631
 632   return Raw2ClustersChamber(rawReader);
 633
 634 }
 635
 636 //_____________________________________________________________________________
 637 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 638 {
 639   //
 640   // Creates clusters from raw data
 641   //
 642
 643   // Create the digits manager
 644   if (!fDigitsManager){
 645     SetBit(knewDM, kTRUE);
 646     fDigitsManager = new AliTRDdigitsManager(kTRUE);
 647     fDigitsManager->CreateArrays();
 648   }
 649
 650   fDigitsManager->SetUseDictionaries(TestBit(kLabels));
 651
 652   // ----- preparing tracklet output -----
 653   if (fReconstructor->IsWritingTracklets()) {
 654     AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 655     if (!trklLoader) {
 656       //AliInfo("Could not get the tracklets data loader, adding it now!");
 657       trklLoader = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
 658       AliRunLoader::Instance()->GetLoader("TRDLoader")->AddDataLoader(trklLoader);
 659     }
 660     AliTreeLoader *trklTreeLoader = dynamic_cast<AliTreeLoader*> (trklLoader->GetBaseLoader("tracklets-raw"));
 661     if (!trklTreeLoader) {
 662       trklTreeLoader = new AliTreeLoader("tracklets-raw", trklLoader);
 663       trklLoader->AddBaseLoader(trklTreeLoader);
 664     }
 665     if (!trklTreeLoader->Tree())
 666       trklTreeLoader->MakeTree();
 667   }
 668
 669   if(!fRawStream)
 670     fRawStream = new AliTRDrawStream(rawReader);
 671   else
 672     fRawStream->SetReader(rawReader);
 673
 674   //if(fReconstructor->IsHLT()){
 675     fRawStream->DisableErrorStorage();
 676   //}
 677
 678   // register tracklet array for output
 679   fRawStream->SetTrackletArray(TrackletsArray("AliTRDtrackletMCM"));
 680   fRawStream->SetTrackArray(TracksArray());
 681
 682   UInt_t det = 0;
 683   while ((det = fRawStream->NextChamber(fDigitsManager)) < AliTRDgeometry::kNdet){
 684     if (fDigitsManager->GetIndexes(det)->HasEntry()){
 685       MakeClusters(det);
 686       fDigitsManager->ClearArrays(det);
 687     }
 688   }
 689
 690   if (fReconstructor->IsWritingTracklets()) {
 691     if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
 692       if (trklLoader) {
 693         if (AliTreeLoader *trklTreeLoader = (AliTreeLoader*) trklLoader->GetBaseLoader("tracklets-raw"))
 694           trklTreeLoader->WriteData("OVERWRITE");
 695         trklLoader->UnloadAll();
 696       }
 697     }
 698   }
 699
 700   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 701
 702   if(!TestBit(knewDM)){
 703     delete fDigitsManager;
 704     fDigitsManager = NULL;
 705     delete fRawStream;
 706     fRawStream = NULL;
 707   }
 708
 709   AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
 710   return kTRUE;
 711
 712 }
 713
 714 //_____________________________________________________________________________
 715 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 716 {
 717   //
 718   // Check if a pad is masked
 719   //
 720
 721   UChar_t status = 0;
 722
 723   if(signal>0 && TESTBIT(signal, 10)){
 724     CLRBIT(signal, 10);
 725     for(int ibit=0; ibit<4; ibit++){
 726       if(TESTBIT(signal, 11+ibit)){
 727         SETBIT(status, ibit);
 728         CLRBIT(signal, 11+ibit);
 729       }
 730     }
 731   }
 732   return status;
 733 }
 734
 735 //_____________________________________________________________________________
 736 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
 737   //
 738   // Set the pad status into out
 739   // First three bits are needed for the position encoding
 740   //
 741   out |= status << 3;
 742 }
 743
 744 //_____________________________________________________________________________
 745 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
 746   //
 747   // return the staus encoding of the corrupted pad
 748   //
 749   return static_cast<UChar_t>(encoding >> 3);
 750 }
 751
 752 //_____________________________________________________________________________
 753 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
 754   //
 755   // Return the position of the corruption
 756   //
 757   return encoding & 7;
 758 }
 759
 760 //_____________________________________________________________________________
 761 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 762 {
 763   //
 764   // Generates the cluster.
 765   //
 766
 767   // Get the digits
 768   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
 769   fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline(det);
 770
 771   // This is to take care of switched off super modules
 772   if (!fDigits->HasData()) return kFALSE;
 773
 774   fIndexes = fDigitsManager->GetIndexes(det);
 775   if (fIndexes->IsAllocated() == kFALSE) {
 776     AliError("Indexes do not exist!");
 777     return kFALSE;
 778   }
 779
 780   AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance();
 781   if (!calibration) {
 782     AliFatal("No AliTRDcalibDB instance available\n");
 783     return kFALSE;
 784   }
 785
 786   if (!fReconstructor){
 787     AliError("Reconstructor not set\n");
 788     return kFALSE;
 789   }
 790
 791   const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
 792
 793   fMaxThresh            = (Short_t)recoParam->GetClusMaxThresh();
 794   fMaxThreshTest        = (Short_t)(recoParam->GetClusMaxThresh()/2+fBaseline);
 795   fSigThresh            = (Short_t)recoParam->GetClusSigThresh();
 796   fMinMaxCutSigma       = recoParam->GetMinMaxCutSigma();
 797   fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma();
 798   const Int_t iEveryNTB = recoParam->GetRecEveryNTB();
 799
 800   Int_t istack  = fIndexes->GetStack();
 801   fLayer  = fIndexes->GetLayer();
 802   Int_t isector = fIndexes->GetSM();
 803
 804   // Start clustering in the chamber
 805
 806   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 807   if (fDet != det) {
 808     AliError("Strange Detector number Missmatch!");
 809     return kFALSE;
 810   }
 811
 812   AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
 813
 814   // TRD space point transformation
 815   fTransform->SetDetector(det);
 816
 817   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 818   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 819   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 820
 821   fColMax    = fDigits->GetNcol();
 822   fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
 823
 824   // Check consistency between OCDB and raw data
 825   Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
 826   if(fReconstructor->IsHLT()){
 827     if((nTimeOCDB > -1) && (fTimeTotal != nTimeOCDB)){
 828       AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value"
 829                       ,fTimeTotal,nTimeOCDB));
 830     }
 831   }else{
 832     if(nTimeOCDB == -1){
 833       AliWarning("Undefined number of timebins in OCDB, using value from raw data.");
 834       if(!fTimeTotal>0){
 835         AliError("Number of timebins in raw data is negative, skipping chamber!");
 836         return kFALSE;
 837       }
 838     }else if(nTimeOCDB == -2){
 839       AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!");
 840       return kFALSE;
 841     }else if(fTimeTotal != nTimeOCDB){
 842       AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber!"
 843                     ,fTimeTotal,nTimeOCDB));
 844       return kFALSE;
 845     }
 846   }
 847
 848   // Detector wise calibration object for the gain factors
 849   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 850   // Calibration object with pad wise values for the gain factors
 851   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 852   // Calibration value for chamber wise gain factor
 853   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 854
 855   // Detector wise calibration object for the noise
 856   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 857   // Calibration object with pad wise values for the noise
 858   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 859   // Calibration value for chamber wise noise
 860   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 861
 862   // Calibration object with the pad status
 863   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
 864   // Calibration object of the online gain
 865   fCalOnGainROC          = calibration->GetOnlineGainTableROC(fDet);
 866
 867   firstClusterROC = -1;
 868   fClusterROC     =  0;
 869
 870   SetBit(kLUT, recoParam->UseLUT());
 871   SetBit(kGAUS, recoParam->UseGAUS());
 872
 873   // Apply the gain and the tail cancelation via digital filter
 874   // Use the configuration from the DCS to find out whether online
 875   // tail cancellation was applied
 876   if(!calibration->HasOnlineTailCancellation()) TailCancelation(recoParam);
 877
 878   MaxStruct curr, last;
 879   Int_t nMaximas = 0, nCorrupted = 0;
 880
 881   // Here the clusterfining is happening
 882
 883   for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
 884     while(fIndexes->NextRCIndex(curr.row, curr.col)){
 885       if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
 886         if(last.row>-1){
 887           if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
 888           CreateCluster(last);
 889         }
 890         last=curr; curr.fivePad=kFALSE;
 891       }
 892     }
 893   }
 894   if(last.row>-1) CreateCluster(last);
 895
 896   if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
 897     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
 898     (*fDebugStream) << "MakeClusters"
 899       << "Detector="   << det
 900       << "NMaxima="    << nMaximas
 901       << "NClusters="  << fClusterROC
 902       << "NCorrupted=" << nCorrupted
 903       << "\n";
 904   }
 905   // if (TestBit(kLabels)) AddLabels();
 906
 907   return kTRUE;
 908
 909 }
 910
 911 //_____________________________________________________________________________
 912 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
 913 {
 914   //
 915   // Returns true if this row,col,time combination is a maximum.
 916   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 917   //
 918
 919   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
 920   Float_t ongain = fCalOnGainROC->GetGainCorrectionFactor(Max.row,Max.col);
 921   Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) * ongain / gain + 0.5f);
 922   if(Signals[1] <= fMaxThresh) return kFALSE;
 923
 924   if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
 925
 926   Short_t noiseMiddleThresh = (Short_t)(fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row));
 927   if (Signals[1] <= noiseMiddleThresh) return kFALSE;
 928
 929   UChar_t status[3]={
 930     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
 931    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
 932    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
 933   };
 934
 935   Short_t signal(0);
 936   if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
 937     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
 938     ongain = fCalOnGainROC->GetGainCorrectionFactor(Max.row,Max.col-1);
 939     Signals[0] = (Short_t)((signal - fBaseline) * ongain / gain + 0.5f);
 940   } else Signals[0] = 0;
 941   if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
 942     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
 943     ongain = fCalOnGainROC->GetGainCorrectionFactor(Max.row,Max.col+1);
 944     Signals[2] = (Short_t)((signal - fBaseline) * ongain / gain + 0.5f);
 945   } else Signals[2] = 0;
 946
 947   if(!(status[0] | status[1] | status[2])) {//all pads are good
 948     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 949       if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
 950         if(Signals[0]<0)Signals[0]=0;
 951         if(Signals[2]<0)Signals[2]=0;
 952         Short_t noiseSumThresh = (Short_t)(fMinLeftRightCutSigma * fCalNoiseDetValue
 953                                            * fCalNoiseROC->GetValue(Max.col, Max.row));
 954         if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
 955         padStatus = 0;
 956         return kTRUE;
 957       }
 958     }
 959   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 960     if(Signals[0]<0)Signals[0]=0;
 961     if(Signals[2]<0)Signals[2]=0;
 962     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
 963       Signals[2]=0;
 964       SetPadStatus(status[2], padStatus);
 965       return kTRUE;
 966     }
 967     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
 968       Signals[0]=0;
 969       SetPadStatus(status[0], padStatus);
 970       return kTRUE;
 971     }
 972     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
 973       Signals[1] = fMaxThresh;
 974       SetPadStatus(status[1], padStatus);
 975       return kTRUE;
 976     }
 977   }
 978   return kFALSE;
 979 }
 980
 981 //_____________________________________________________________________________
 982 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
 983 {
 984   //
 985   // Look for 5 pad cluster with minimum in the middle
 986   // Gives back the ratio
 987   //
 988
 989   if (ThisMax.col >= fColMax - 3) return kFALSE;
 990   Float_t gain;
 991   if (ThisMax.col < fColMax - 5){
 992     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
 993     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
 994       return kFALSE;
 995   }
 996   if (ThisMax.col > 1) {
 997     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
 998     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
 999       return kFALSE;
1000   }
1001
1002   const Float_t kEpsilon = 0.01;
1003   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
1004       NeighbourMax.signals[1], NeighbourMax.signals[2]};
1005
1006   // Unfold the two maxima and set the signal on
1007   // the overlapping pad to the ratio
1008   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
1009   ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
1010   NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
1011   ThisMax.fivePad=kTRUE;
1012   NeighbourMax.fivePad=kTRUE;
1013   return kTRUE;
1014
1015 }
1016
1017 //_____________________________________________________________________________
1018 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
1019 {
1020   //
1021   // Creates a cluster at the given position and saves it in fRecPoints
1022   //
1023
1024   Int_t nPadCount = 1;
1025   Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
1026   if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
1027
1028   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
1029   cluster.SetNPads(nPadCount);
1030   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1031   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1032   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
1033
1034   cluster.SetFivePad(Max.fivePad);
1035   // set pads status for the cluster
1036   UChar_t maskPosition = GetCorruption(Max.padStatus);
1037   if (maskPosition) {
1038     cluster.SetPadMaskedPosition(maskPosition);
1039     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1040   }
1041   cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
1042
1043   // Transform the local cluster coordinates into calibrated
1044   // space point positions defined in the local tracking system.
1045   // Here the calibration for T0, Vdrift and ExB is applied as well.
1046   if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
1047
1048   // Temporarily store the Max.Row, column and time bin of the center pad
1049   // Used to later on assign the track indices
1050   cluster.SetLabel(Max.row, 0);
1051   cluster.SetLabel(Max.col, 1);
1052   cluster.SetLabel(Max.time,2);
1053
1054   //needed for HLT reconstruction
1055   AddClusterToArray(&cluster);
1056
1057   // Store the index of the first cluster in the current ROC
1058   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1059
1060   fNoOfClusters++;
1061   fClusterROC++;
1062 }
1063
1064 //_____________________________________________________________________________
1065 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1066 {
1067 // Calculate number of pads/cluster and
1068 // ADC signals at position 0, 1, 5 and 6
1069
1070   Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
1071   Float_t gain(1.); Short_t signal(0);
1072   // Store the amplitudes of the pads in the cluster for later analysis
1073   // and check whether one of these pads is masked in the database
1074   signals[3]=Max.signals[1];
1075   Int_t ipad(1), jpad(0);
1076   // Look to the right
1077   while((jpad = Max.col-ipad)){
1078     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1079     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1080     tmp = (signal - fBaseline) / gain + 0.5f;
1081     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1082     if(signal<fSigThresh) break; // signal under threshold
1083     nPadCount++;
1084     if(ipad<=3) signals[3 - ipad] = signal;
1085     ipad++;
1086   }
1087   ipad=1;
1088   // Look to the left
1089   while((jpad = Max.col+ipad)<fColMax){
1090     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1091     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1092     tmp = (signal - fBaseline) / gain + 0.5f;
1093     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1094     if(signal<fSigThresh) break; // signal under threshold
1095     nPadCount++;
1096     if(ipad<=3) signals[3 + ipad] = signal;
1097     ipad++;
1098   }
1099
1100   AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1101     , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
1102 }
1103
1104 //_____________________________________________________________________________
1105 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1106 {
1107   //
1108   // Add a cluster to the array
1109   //
1110
1111   Int_t n = RecPoints()->GetEntriesFast();
1112   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1113   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1114 }
1115
1116 //_____________________________________________________________________________
1117 Bool_t AliTRDclusterizer::AddLabels()
1118 {
1119   //
1120   // Add the track indices to the found clusters
1121   //
1122
1123   const Int_t   kNclus  = 3;
1124   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1125   const Int_t   kNtrack = kNdict * kNclus;
1126
1127   Int_t iClusterROC = 0;
1128
1129   Int_t row  = 0;
1130   Int_t col  = 0;
1131   Int_t time = 0;
1132   Int_t iPad = 0;
1133
1134   // Temporary array to collect the track indices
1135   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1136
1137   // Loop through the dictionary arrays one-by-one
1138   // to keep memory consumption low
1139   AliTRDarrayDictionary *tracksIn = 0;  //mod
1140   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1141
1142     // tracksIn should be expanded beforehand!
1143     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1144
1145     // Loop though the clusters found in this ROC
1146     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1147
1148       AliTRDcluster *cluster = (AliTRDcluster *)
1149                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1150       row  = cluster->GetLabel(0);
1151       col  = cluster->GetLabel(1);
1152       time = cluster->GetLabel(2);
1153
1154       for (iPad = 0; iPad < kNclus; iPad++) {
1155         Int_t iPadCol = col - 1 + iPad;
1156         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1157         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1158       }
1159
1160     }
1161
1162   }
1163
1164   // Copy the track indices into the cluster
1165   // Loop though the clusters found in this ROC
1166   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1167
1168     AliTRDcluster *cluster = (AliTRDcluster *)
1169       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1170     cluster->SetLabel(-9999,0);
1171     cluster->SetLabel(-9999,1);
1172     cluster->SetLabel(-9999,2);
1173
1174     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1175
1176   }
1177
1178   delete [] idxTracks;
1179
1180   return kTRUE;
1181
1182 }
1183
1184 //_____________________________________________________________________________
1185 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1186 {
1187   //
1188   // Method to unfold neighbouring maxima.
1189   // The charge ratio on the overlapping pad is calculated
1190   // until there is no more change within the range given by eps.
1191   // The resulting ratio is then returned to the calling method.
1192   //
1193
1194   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1195   if (!calibration) {
1196     AliError("No AliTRDcalibDB instance available\n");
1197     return kFALSE;
1198   }
1199
1200   Int_t   irc                = 0;
1201   Int_t   itStep             = 0;                 // Count iteration steps
1202
1203   Double_t ratio             = 0.5;               // Start value for ratio
1204   Double_t prevRatio         = 0.0;               // Store previous ratio
1205
1206   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1207   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1208   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1209
1210   // Start the iteration
1211   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1212
1213     itStep++;
1214     prevRatio = ratio;
1215
1216     // Cluster position according to charge ratio
1217     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1218                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1219     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1220                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1221
1222     // Set cluster charge ratio
1223     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1224     Double_t ampLeft  = padSignal[1] / newSignal[1];
1225     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1226     Double_t ampRight = padSignal[3] / newSignal[1];
1227
1228     // Apply pad response to parameters
1229     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1230     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1231
1232     // Calculate new overlapping ratio
1233     ratio = TMath::Min((Double_t) 1.0
1234                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1235
1236   }
1237
1238   return ratio;
1239
1240 }
1241
1242 //_____________________________________________________________________________
1243 void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1244 {
1245   //
1246   // Applies the tail cancelation
1247   //
1248
1249   Int_t nexp = recoParam->GetTCnexp();
1250   if(!nexp) return;
1251
1252   Int_t iRow  = 0;
1253   Int_t iCol  = 0;
1254   Int_t iTime = 0;
1255
1256   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1257   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1258   while(fIndexes->NextRCIndex(iRow, iCol))
1259     {
1260       // if corrupted then don't make the tail cancallation
1261       if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1262
1263       if(debugStreaming){
1264         for (iTime = 0; iTime < fTimeTotal; iTime++)
1265           (*fDebugStream) << "TailCancellation"
1266                           << "col="  << iCol
1267                           << "row="  << iRow
1268                           << "\n";
1269       }
1270
1271       // Apply the tail cancelation via the digital filter
1272       //DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1273       ApplyTCTM(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1274     } // while irow icol
1275
1276   return;
1277
1278 }
1279
1280
1281 //_____________________________________________________________________________
1282 void AliTRDclusterizer::ApplyTCTM(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1283 {
1284   //
1285   // Steer tail cancellation
1286   //
1287
1288
1289   switch(nexp) {
1290   case 1:
1291   case 2:
1292     DeConvExp(arr,nTime,nexp);
1293     break;
1294   case -1:
1295     ConvExp(arr,nTime);
1296     break;
1297   case -2:
1298     DeConvExp(arr,nTime,1);
1299     ConvExp(arr,nTime);
1300     break;
1301   default:
1302     break;
1303   }
1304 }
1305
1306
1307 //_____________________________________________________________________________
1308 void AliTRDclusterizer::ConvExp(Short_t *const arr, const Int_t nTime)
1309 {
1310   //
1311   // Tail maker
1312   //
1313
1314   // Initialization (coefficient = alpha, rates = lambda)
1315   Float_t slope = 1.0;
1316   Float_t coeff = 0.5;
1317   Float_t rate;
1318
1319   Double_t par[4];
1320   fReconstructor->GetRecoParam()->GetTCParams(par);
1321   slope = par[1];
1322   coeff = par[3];
1323
1324   Double_t dt = 0.1;
1325
1326   rate = TMath::Exp(-dt/(slope));
1327
1328   Float_t reminder =  .0;
1329   Float_t correction = 0.0;
1330   Float_t result     = 0.0;
1331
1332   for (int i = nTime-1; i >= 0; i--) {
1333
1334     result = arr[i] + correction - fBaseline;    // No rescaling
1335     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1336
1337     correction = 0.0;
1338
1339     correction += reminder = rate * (reminder + coeff * result);
1340   }
1341 }
1342
1343
1344 //_____________________________________________________________________________
1345 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1346 {
1347   //
1348   // Tail cancellation by deconvolution for PASA v4 TRF
1349   //
1350
1351   Float_t rates[2];
1352   Float_t coefficients[2];
1353
1354   // Initialization (coefficient = alpha, rates = lambda)
1355   Float_t r1 = 1.0;
1356   Float_t r2 = 1.0;
1357   Float_t c1 = 0.5;
1358   Float_t c2 = 0.5;
1359
1360   if (nexp == 1) {   // 1 Exponentials
1361     r1 = 1.156;
1362     r2 = 0.130;
1363     c1 = 0.066;
1364     c2 = 0.000;
1365   }
1366   if (nexp == 2) {   // 2 Exponentials
1367     Double_t par[4];
1368     fReconstructor->GetRecoParam()->GetTCParams(par);
1369     r1 = par[0];//1.156;
1370     r2 = par[1];//0.130;
1371     c1 = par[2];//0.114;
1372     c2 = par[3];//0.624;
1373   }
1374
1375   coefficients[0] = c1;
1376   coefficients[1] = c2;
1377
1378   Double_t dt = 0.1;
1379
1380   rates[0] = TMath::Exp(-dt/(r1));
1381   rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
1382
1383   Float_t reminder[2] = { .0, .0 };
1384   Float_t correction = 0.0;
1385   Float_t result     = 0.0;
1386
1387   for (int i = 0; i < nTime; i++) {
1388
1389     result = arr[i] - correction - fBaseline;    // No rescaling
1390     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1391
1392     correction = 0.0;
1393     for (int k = 0; k < 2; k++) {
1394       correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1395     }
1396
1397   }
1398
1399 }
1400
1401 //_____________________________________________________________________________
1402 void AliTRDclusterizer::ResetRecPoints()
1403 {
1404   //
1405   // Resets the list of rec points
1406   //
1407
1408   if (fRecPoints) {
1409     fRecPoints->Clear();
1410     fNoOfClusters = 0;
1411     //    delete fRecPoints;
1412   }
1413 }
1414
1415 //_____________________________________________________________________________
1416 TClonesArray *AliTRDclusterizer::RecPoints()
1417 {
1418   //
1419   // Returns the list of rec points
1420   //
1421
1422   if (!fRecPoints) {
1423     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1424       // determine number of clusters which has to be allocated
1425       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1426
1427       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1428     }
1429     //SetClustersOwner(kTRUE);
1430     AliTRDReconstructor::SetClusters(0x0);
1431   }
1432   return fRecPoints;
1433
1434 }
1435
1436 //_____________________________________________________________________________
1437 TClonesArray *AliTRDclusterizer::TrackletsArray(const TString &trkltype)
1438 {
1439   //
1440   // Returns the array of on-line tracklets
1441   //
1442
1443   if (trkltype.Length() != 0) {
1444     if (!fTracklets) {
1445       fTracklets = new TClonesArray(trkltype, 200);
1446   }
1447     else if (TClass::GetClass(trkltype.Data()) != fTracklets->GetClass()){
1448       fTracklets->Delete();
1449       delete fTracklets;
1450       fTracklets = new TClonesArray(trkltype, 200);
1451     }
1452   }
1453   return fTracklets;
1454 }
1455
1456 //_____________________________________________________________________________
1457 TClonesArray* AliTRDclusterizer::TracksArray()
1458 {
1459   // return array of GTU tracks (create TClonesArray if necessary)
1460
1461   if (!fTracks) {
1462     fTracks = new TClonesArray("AliESDTrdTrack",100);
1463   }
1464   return fTracks;
1465 }