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