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