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