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   ,fCalOnGainROC(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   ,fCalOnGainROC(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   ,fCalOnGainROC(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   fCalOnGainROC          = calibration->GetOnlineGainTableROC(fDet);
 901
 902   firstClusterROC = -1;
 903   fClusterROC     =  0;
 904
 905   SetBit(kLUT, recoParam->UseLUT());
 906   SetBit(kGAUS, recoParam->UseGAUS());
 907
 908   // Apply the gain and the tail cancelation via digital filter
 909   // Use the configuration from the DCS to find out whether online
 910   // tail cancellation was applied
 911   if(!calibration->HasOnlineTailCancellation()){
 912     // save a copy of raw data
 913     if(TestBit(kRawSignal)){
 914       if(fDigitsRaw){
 915         fDigitsRaw->~AliTRDarrayADC();
 916         new(fDigitsRaw) AliTRDarrayADC(*fDigits);
 917       } else fDigitsRaw = new AliTRDarrayADC(*fDigits);
 918     }
 919     TailCancelation(recoParam);
 920   }
 921
 922   MaxStruct curr, last;
 923   Int_t nMaximas = 0, nCorrupted = 0;
 924
 925   // Here the clusterfining is happening
 926
 927   for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
 928     while(fIndexes->NextRCIndex(curr.row, curr.col)){
 929       if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
 930         if(last.row>-1){
 931           if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
 932           CreateCluster(last);
 933         }
 934         last=curr; curr.fivePad=kFALSE;
 935       }
 936     }
 937   }
 938   if(last.row>-1) CreateCluster(last);
 939
 940   if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
 941     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
 942     (*fDebugStream) << "MakeClusters"
 943       << "Detector="   << det
 944       << "NMaxima="    << nMaximas
 945       << "NClusters="  << fClusterROC
 946       << "NCorrupted=" << nCorrupted
 947       << "\n";
 948   }
 949   if (TestBit(kLabels)) AddLabels();
 950
 951   return kTRUE;
 952
 953 }
 954
 955 //_____________________________________________________________________________
 956 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Float_t *const Signals)
 957 {
 958   //
 959   // Returns true if this row,col,time combination is a maximum.
 960   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 961   //
 962
 963   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
 964   Float_t ongain = fCalOnGainROC ? fCalOnGainROC->GetGainCorrectionFactor(Max.row,Max.col) : 1;
 965   Signals[1] = (fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) * ongain / gain + 0.5f;
 966   if(Signals[1] <= fMaxThresh) return kFALSE;
 967
 968   if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
 969
 970   Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
 971   if (Signals[1] <= noiseMiddleThresh) return kFALSE;
 972
 973   UChar_t status[3]={
 974     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
 975    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
 976    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
 977   };
 978
 979   Short_t signal(0);
 980   if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
 981     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
 982     ongain = fCalOnGainROC ? fCalOnGainROC->GetGainCorrectionFactor(Max.row,Max.col-1) : 1;
 983     Signals[0] = (signal - fBaseline) * ongain / gain + 0.5f;
 984   } else Signals[0] = 0.;
 985   if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
 986     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
 987     ongain = fCalOnGainROC ? fCalOnGainROC->GetGainCorrectionFactor(Max.row,Max.col+1) : 1;
 988     Signals[2] = (signal - fBaseline) * ongain / gain + 0.5f;
 989   } else Signals[2] = 0.;
 990
 991   if(!(status[0] | status[1] | status[2])) {//all pads are good
 992     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 993       if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
 994         if(Signals[0]<0) Signals[0]=0.;
 995         if(Signals[2]<0) Signals[2]=0.;
 996         Float_t noiseSumThresh = fMinLeftRightCutSigma * fCalNoiseDetValue
 997                                * fCalNoiseROC->GetValue(Max.col, Max.row);
 998         if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
 999         padStatus = 0;
1000         return kTRUE;
1001       }
1002     }
1003   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
1004     if(Signals[0]<0)Signals[0]=0;
1005     if(Signals[2]<0)Signals[2]=0;
1006     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
1007       Signals[2]=0;
1008       SetPadStatus(status[2], padStatus);
1009       return kTRUE;
1010     }
1011     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
1012       Signals[0]=0;
1013       SetPadStatus(status[0], padStatus);
1014       return kTRUE;
1015     }
1016     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
1017       Signals[1] = fMaxThresh;
1018       SetPadStatus(status[1], padStatus);
1019       return kTRUE;
1020     }
1021   }
1022   return kFALSE;
1023 }
1024
1025 //_____________________________________________________________________________
1026 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
1027 {
1028   //
1029   // Look for 5 pad cluster with minimum in the middle
1030   // Gives back the ratio
1031   //
1032
1033   if (ThisMax.col >= fColMax - 3) return kFALSE;
1034   Float_t gain;
1035   if (ThisMax.col < fColMax - 5){
1036     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
1037     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
1038       return kFALSE;
1039   }
1040   if (ThisMax.col > 1) {
1041     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
1042     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
1043       return kFALSE;
1044   }
1045
1046   const Float_t kEpsilon = 0.01;
1047   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
1048       NeighbourMax.signals[1], NeighbourMax.signals[2]};
1049
1050   // Unfold the two maxima and set the signal on
1051   // the overlapping pad to the ratio
1052   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
1053   ThisMax.signals[2] = ThisMax.signals[2]*ratio + 0.5f;
1054   NeighbourMax.signals[0] = NeighbourMax.signals[0]*(1-ratio) + 0.5f;
1055   ThisMax.fivePad=kTRUE;
1056   NeighbourMax.fivePad=kTRUE;
1057   return kTRUE;
1058
1059 }
1060
1061 //_____________________________________________________________________________
1062 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
1063 {
1064   //
1065   // Creates a cluster at the given position and saves it in fRecPoints
1066   //
1067
1068   Int_t nPadCount = 1;
1069   Short_t signals[7] = { 0, 0, (Short_t)Max.signals[0], (Short_t)Max.signals[1], (Short_t)Max.signals[2], 0, 0 };
1070   if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
1071
1072   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
1073   cluster.SetNPads(nPadCount);
1074   cluster.SetQ(Max.signals[0]+Max.signals[1]+Max.signals[2]);
1075   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1076   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1077   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
1078
1079   cluster.SetFivePad(Max.fivePad);
1080   // set pads status for the cluster
1081   UChar_t maskPosition = GetCorruption(Max.padStatus);
1082   if (maskPosition) {
1083     cluster.SetPadMaskedPosition(maskPosition);
1084     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1085   }
1086   cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
1087
1088   // Transform the local cluster coordinates into calibrated
1089   // space point positions defined in the local tracking system.
1090   // Here the calibration for T0, Vdrift and ExB is applied as well.
1091   if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
1092   // Store raw signals in cluster. This MUST be called after position reconstruction !
1093   // Xianguo Lu and Alex Bercuci 19.03.2012
1094   if(TestBit(kRawSignal) && fDigitsRaw){
1095     Short_t rawSignal[7] = {0, 0, 0, fDigitsRaw->GetData(Max.row, Max.col, Max.time), 0, 0, 0};
1096     for(Int_t ipad(1); ipad<=3; ipad++){
1097       rawSignal[3 - ipad] = (Max.col-ipad>=0)?fDigitsRaw->GetData(Max.row, Max.col-ipad, Max.time):0; // Look to the left
1098       rawSignal[3 + ipad] = (Max.col+ipad<fColMax)?fDigitsRaw->GetData(Max.row, Max.col+ipad, Max.time):0; // Look to the right
1099     }
1100     cluster.SetSignals(rawSignal, kTRUE);
1101   }
1102   // Temporarily store the Max.Row, column and time bin of the center pad
1103   // Used to later on assign the track indices
1104   cluster.SetLabel(Max.row, 0);
1105   cluster.SetLabel(Max.col, 1);
1106   cluster.SetLabel(Max.time,2);
1107
1108   //needed for HLT reconstruction
1109   AddClusterToArray(&cluster);
1110
1111   // Store the index of the first cluster in the current ROC
1112   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1113
1114   fNoOfClusters++;
1115   fClusterROC++;
1116 }
1117
1118 //_____________________________________________________________________________
1119 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1120 {
1121 // Calculate number of pads/cluster and
1122 // ADC signals at position 0, 1, 5 and 6
1123
1124   Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
1125   Float_t gain(1.); Short_t signal(0);
1126   // Store the amplitudes of the pads in the cluster for later analysis
1127   // and check whether one of these pads is masked in the database
1128   signals[3]=Max.signals[1];
1129   Int_t ipad(1), jpad(0);
1130   // Look to the right
1131   while((jpad = Max.col-ipad)){
1132     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1133     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1134     tmp = (signal - fBaseline) / gain + 0.5f;
1135     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1136     if(signal<fSigThresh) break; // signal under threshold
1137     nPadCount++;
1138     if(ipad<=3) signals[3 - ipad] = signal;
1139     ipad++;
1140   }
1141   ipad=1;
1142   // Look to the left
1143   while((jpad = Max.col+ipad)<fColMax){
1144     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1145     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1146     tmp = (signal - fBaseline) / gain + 0.5f;
1147     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1148     if(signal<fSigThresh) break; // signal under threshold
1149     nPadCount++;
1150     if(ipad<=3) signals[3 + ipad] = signal;
1151     ipad++;
1152   }
1153
1154   AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1155     , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
1156 }
1157
1158 //_____________________________________________________________________________
1159 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1160 {
1161   //
1162   // Add a cluster to the array
1163   //
1164
1165   Int_t n = RecPoints()->GetEntriesFast();
1166   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1167   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1168 }
1169
1170 //_____________________________________________________________________________
1171 Bool_t AliTRDclusterizer::AddLabels()
1172 {
1173   //
1174   // Add the track indices to the found clusters
1175   //
1176
1177   const Int_t   kNclus  = 3;
1178   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1179   const Int_t   kNtrack = kNdict * kNclus;
1180
1181   Int_t iClusterROC = 0;
1182
1183   Int_t row  = 0;
1184   Int_t col  = 0;
1185   Int_t time = 0;
1186   Int_t iPad = 0;
1187
1188   // Temporary array to collect the track indices
1189   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1190
1191   // Loop through the dictionary arrays one-by-one
1192   // to keep memory consumption low
1193   AliTRDarrayDictionary *tracksIn(NULL);  //mod
1194   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1195
1196     // tracksIn should be expanded beforehand!
1197     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1198
1199     // Loop though the clusters found in this ROC
1200     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1201
1202       AliTRDcluster *cluster = (AliTRDcluster *)
1203                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1204       row  = cluster->GetLabel(0);
1205       col  = cluster->GetLabel(1);
1206       time = cluster->GetLabel(2);
1207
1208       for (iPad = 0; iPad < kNclus; iPad++) {
1209         Int_t iPadCol = col - 1 + iPad;
1210         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1211         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1212       }
1213
1214     }
1215
1216   }
1217
1218   // Copy the track indices into the cluster
1219   // Loop though the clusters found in this ROC
1220   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1221
1222     AliTRDcluster *cluster = (AliTRDcluster *)
1223       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1224     cluster->SetLabel(-9999,0);
1225     cluster->SetLabel(-9999,1);
1226     cluster->SetLabel(-9999,2);
1227
1228     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1229
1230   }
1231
1232   delete [] idxTracks;
1233
1234   return kTRUE;
1235
1236 }
1237
1238 //_____________________________________________________________________________
1239 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1240 {
1241   //
1242   // Method to unfold neighbouring maxima.
1243   // The charge ratio on the overlapping pad is calculated
1244   // until there is no more change within the range given by eps.
1245   // The resulting ratio is then returned to the calling method.
1246   //
1247
1248   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1249   if (!calibration) {
1250     AliError("No AliTRDcalibDB instance available\n");
1251     return kFALSE;
1252   }
1253
1254   Int_t   irc                = 0;
1255   Int_t   itStep             = 0;                 // Count iteration steps
1256
1257   Double_t ratio             = 0.5;               // Start value for ratio
1258   Double_t prevRatio         = 0.0;               // Store previous ratio
1259
1260   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1261   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1262   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1263
1264   // Start the iteration
1265   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1266
1267     itStep++;
1268     prevRatio = ratio;
1269
1270     // Cluster position according to charge ratio
1271     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1272                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1273     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1274                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1275
1276     // Set cluster charge ratio
1277     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1278     Double_t ampLeft  = padSignal[1] / newSignal[1];
1279     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1280     Double_t ampRight = padSignal[3] / newSignal[1];
1281
1282     // Apply pad response to parameters
1283     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1284     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1285
1286     // Calculate new overlapping ratio
1287     ratio = TMath::Min((Double_t) 1.0
1288                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1289
1290   }
1291
1292   return ratio;
1293
1294 }
1295
1296 //_____________________________________________________________________________
1297 void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1298 {
1299   //
1300   // Applies the tail cancelation
1301   //
1302
1303   Int_t nexp = recoParam->GetTCnexp();
1304   if(!nexp) return;
1305
1306   Int_t iRow  = 0;
1307   Int_t iCol  = 0;
1308   Int_t iTime = 0;
1309
1310   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1311   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1312   while(fIndexes->NextRCIndex(iRow, iCol))
1313     {
1314       // if corrupted then don't make the tail cancallation
1315       if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1316
1317       if(debugStreaming){
1318         for (iTime = 0; iTime < fTimeTotal; iTime++)
1319           (*fDebugStream) << "TailCancellation"
1320                           << "col="  << iCol
1321                           << "row="  << iRow
1322                           << "\n";
1323       }
1324
1325       // Apply the tail cancelation via the digital filter
1326       //DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1327       ApplyTCTM(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1328     } // while irow icol
1329
1330   return;
1331
1332 }
1333
1334
1335 //_____________________________________________________________________________
1336 void AliTRDclusterizer::ApplyTCTM(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1337 {
1338   //
1339   // Steer tail cancellation
1340   //
1341
1342
1343   switch(nexp) {
1344   case 1:
1345   case 2:
1346     DeConvExp(arr,nTime,nexp);
1347     break;
1348   case -1:
1349     ConvExp(arr,nTime);
1350     break;
1351   case -2:
1352     DeConvExp(arr,nTime,1);
1353     ConvExp(arr,nTime);
1354     break;
1355   default:
1356     break;
1357   }
1358 }
1359
1360
1361 //_____________________________________________________________________________
1362 void AliTRDclusterizer::ConvExp(Short_t *const arr, const Int_t nTime)
1363 {
1364   //
1365   // Tail maker
1366   //
1367
1368   // Initialization (coefficient = alpha, rates = lambda)
1369   Float_t slope = 1.0;
1370   Float_t coeff = 0.5;
1371   Float_t rate;
1372
1373   Double_t par[4];
1374   fReconstructor->GetRecoParam()->GetTCParams(par);
1375   slope = par[1];
1376   coeff = par[3];
1377
1378   Double_t dt = 0.1;
1379
1380   rate = TMath::Exp(-dt/(slope));
1381
1382   Float_t reminder =  .0;
1383   Float_t correction = 0.0;
1384   Float_t result     = 0.0;
1385
1386   for (int i = nTime-1; i >= 0; i--) {
1387
1388     result = arr[i] + correction - fBaseline;    // No rescaling
1389     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1390
1391     correction = 0.0;
1392
1393     correction += reminder = rate * (reminder + coeff * result);
1394   }
1395 }
1396
1397
1398 //_____________________________________________________________________________
1399 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1400 {
1401   //
1402   // Tail cancellation by deconvolution for PASA v4 TRF
1403   //
1404
1405   Float_t rates[2];
1406   Float_t coefficients[2];
1407
1408   // Initialization (coefficient = alpha, rates = lambda)
1409   Float_t r1 = 1.0;
1410   Float_t r2 = 1.0;
1411   Float_t c1 = 0.5;
1412   Float_t c2 = 0.5;
1413
1414   if (nexp == 1) {   // 1 Exponentials
1415     r1 = 1.156;
1416     r2 = 0.130;
1417     c1 = 0.066;
1418     c2 = 0.000;
1419   }
1420   if (nexp == 2) {   // 2 Exponentials
1421     Double_t par[4];
1422     fReconstructor->GetRecoParam()->GetTCParams(par);
1423     r1 = par[0];//1.156;
1424     r2 = par[1];//0.130;
1425     c1 = par[2];//0.114;
1426     c2 = par[3];//0.624;
1427   }
1428
1429   coefficients[0] = c1;
1430   coefficients[1] = c2;
1431
1432   Double_t dt = 0.1;
1433
1434   rates[0] = TMath::Exp(-dt/(r1));
1435   rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
1436
1437   Float_t reminder[2] = { .0, .0 };
1438   Float_t correction = 0.0;
1439   Float_t result     = 0.0;
1440
1441   for (int i = 0; i < nTime; i++) {
1442
1443     result = arr[i] - correction - fBaseline;    // No rescaling
1444     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1445
1446     correction = 0.0;
1447     for (int k = 0; k < 2; k++) {
1448       correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1449     }
1450
1451   }
1452
1453 }
1454
1455 //_____________________________________________________________________________
1456 void AliTRDclusterizer::ResetRecPoints()
1457 {
1458   //
1459   // Resets the list of rec points
1460   //
1461
1462   if (fRecPoints) {
1463     fRecPoints->Clear();
1464     fNoOfClusters = 0;
1465     //    delete fRecPoints;
1466   }
1467 }
1468
1469 //_____________________________________________________________________________
1470 TClonesArray *AliTRDclusterizer::RecPoints()
1471 {
1472   //
1473   // Returns the list of rec points
1474   //
1475
1476   if (!fRecPoints) {
1477     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1478       // determine number of clusters which has to be allocated
1479       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1480
1481       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1482     }
1483     //SetClustersOwner(kTRUE);
1484     AliTRDReconstructor::SetClusters(0x0);
1485   }
1486   return fRecPoints;
1487
1488 }
1489
1490 //_____________________________________________________________________________
1491 TClonesArray *AliTRDclusterizer::TrackletsArray(const TString &trkltype)
1492 {
1493   //
1494   // Returns the array of on-line tracklets
1495   //
1496
1497   if (trkltype.Length() != 0) {
1498     if (!fTracklets) {
1499       fTracklets = new TClonesArray(trkltype, 200);
1500   }
1501     else if (TClass::GetClass(trkltype.Data()) != fTracklets->GetClass()){
1502       fTracklets->Delete();
1503       delete fTracklets;
1504       fTracklets = new TClonesArray(trkltype, 200);
1505     }
1506   }
1507   return fTracklets;
1508 }
1509
1510 //_____________________________________________________________________________
1511 TClonesArray* AliTRDclusterizer::TracksArray()
1512 {
1513   // return array of GTU tracks (create TClonesArray if necessary)
1514
1515   if (!fTracks) {
1516     fTracks = new TClonesArray("AliESDTrdTrack",100);
1517   }
1518   return fTracks;
1519 }