TRD/AliTRDclusterizer.cxx

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