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