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