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