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