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(!TestBit(kHLT)) if(!fTransform->Transform(&cluster)) return;
1026
1027   // Temporarily store the Max.Row, column and time bin of the center pad
1028   // Used to later on assign the track indices
1029   cluster.SetLabel(Max.row, 0);
1030   cluster.SetLabel(Max.col, 1);
1031   cluster.SetLabel(Max.time,2);
1032
1033   //needed for HLT reconstruction
1034   AddClusterToArray(&cluster);
1035
1036   // Store the index of the first cluster in the current ROC
1037   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1038
1039   fNoOfClusters++;
1040   fClusterROC++;
1041 }
1042
1043 //_____________________________________________________________________________
1044 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1045 {
1046   // Look to the right
1047   Int_t ii = 1;
1048   while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
1049     nPadCount++;
1050     ii++;
1051     if (Max.col < ii) break;
1052   }
1053   // Look to the left
1054   ii = 1;
1055   while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
1056     nPadCount++;
1057     ii++;
1058     if (Max.col+ii >= fColMax) break;
1059   }
1060
1061   // Store the amplitudes of the pads in the cluster for later analysis
1062   // and check whether one of these pads is masked in the database
1063   signals[2]=Max.signals[0];
1064   signals[3]=Max.signals[1];
1065   signals[4]=Max.signals[2];
1066   Float_t gain;
1067   for(Int_t i = 0; i<2; i++)
1068     {
1069       if(Max.col+i >= 3){
1070         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1071         signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1072       }
1073       if(Max.col+3-i < fColMax){
1074         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1075         signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1076       }
1077     }
1078   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1079     if ((jPad >= 0) && (jPad < fColMax))
1080       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1081       }*/
1082 }
1083
1084 //_____________________________________________________________________________
1085 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1086 {
1087   //
1088   // Add a cluster to the array
1089   //
1090
1091   Int_t n = RecPoints()->GetEntriesFast();
1092   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1093   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1094 }
1095
1096 //_____________________________________________________________________________
1097 void AliTRDclusterizer::AddTrackletsToArray()
1098 {
1099   //
1100   // Add the online tracklets of this chamber to the array
1101   //
1102
1103   UInt_t* trackletword;
1104   for(Int_t side=0; side<2; side++)
1105     {
1106       Int_t trkl=0;
1107       trackletword=fTrackletContainer[side];
1108       while(trackletword[trkl]>0){
1109         Int_t n = TrackletsArray()->GetEntriesFast();
1110         AliTRDtrackletWord tmp(trackletword[trkl]);
1111         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1112         trkl++;
1113       }
1114     }
1115 }
1116
1117 //_____________________________________________________________________________
1118 Bool_t AliTRDclusterizer::AddLabels()
1119 {
1120   //
1121   // Add the track indices to the found clusters
1122   //
1123
1124   const Int_t   kNclus  = 3;
1125   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1126   const Int_t   kNtrack = kNdict * kNclus;
1127
1128   Int_t iClusterROC = 0;
1129
1130   Int_t row  = 0;
1131   Int_t col  = 0;
1132   Int_t time = 0;
1133   Int_t iPad = 0;
1134
1135   // Temporary array to collect the track indices
1136   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1137
1138   // Loop through the dictionary arrays one-by-one
1139   // to keep memory consumption low
1140   AliTRDarrayDictionary *tracksIn = 0;  //mod
1141   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1142
1143     // tracksIn should be expanded beforehand!
1144     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1145
1146     // Loop though the clusters found in this ROC
1147     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1148
1149       AliTRDcluster *cluster = (AliTRDcluster *)
1150                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1151       row  = cluster->GetLabel(0);
1152       col  = cluster->GetLabel(1);
1153       time = cluster->GetLabel(2);
1154
1155       for (iPad = 0; iPad < kNclus; iPad++) {
1156         Int_t iPadCol = col - 1 + iPad;
1157         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1158         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1159       }
1160
1161     }
1162
1163   }
1164
1165   // Copy the track indices into the cluster
1166   // Loop though the clusters found in this ROC
1167   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1168
1169     AliTRDcluster *cluster = (AliTRDcluster *)
1170       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1171     cluster->SetLabel(-9999,0);
1172     cluster->SetLabel(-9999,1);
1173     cluster->SetLabel(-9999,2);
1174
1175     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1176
1177   }
1178
1179   delete [] idxTracks;
1180
1181   return kTRUE;
1182
1183 }
1184
1185 //_____________________________________________________________________________
1186 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1187 {
1188   //
1189   // Method to unfold neighbouring maxima.
1190   // The charge ratio on the overlapping pad is calculated
1191   // until there is no more change within the range given by eps.
1192   // The resulting ratio is then returned to the calling method.
1193   //
1194
1195   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1196   if (!calibration) {
1197     AliError("No AliTRDcalibDB instance available\n");
1198     return kFALSE;
1199   }
1200
1201   Int_t   irc                = 0;
1202   Int_t   itStep             = 0;                 // Count iteration steps
1203
1204   Double_t ratio             = 0.5;               // Start value for ratio
1205   Double_t prevRatio         = 0.0;               // Store previous ratio
1206
1207   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1208   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1209   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1210
1211   // Start the iteration
1212   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1213
1214     itStep++;
1215     prevRatio = ratio;
1216
1217     // Cluster position according to charge ratio
1218     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1219                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1220     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1221                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1222
1223     // Set cluster charge ratio
1224     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1225     Double_t ampLeft  = padSignal[1] / newSignal[1];
1226     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1227     Double_t ampRight = padSignal[3] / newSignal[1];
1228
1229     // Apply pad response to parameters
1230     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1231     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1232
1233     // Calculate new overlapping ratio
1234     ratio = TMath::Min((Double_t) 1.0
1235                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1236
1237   }
1238
1239   return ratio;
1240
1241 }
1242
1243 //_____________________________________________________________________________
1244 void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1245 {
1246   //
1247   // Applies the tail cancelation
1248   //
1249
1250   Int_t iRow  = 0;
1251   Int_t iCol  = 0;
1252   Int_t iTime = 0;
1253
1254   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1255   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1256   Int_t nexp = recoParam->GetTCnexp();
1257   while(fIndexes->NextRCIndex(iRow, iCol))
1258     {
1259       // if corrupted then don't make the tail cancallation
1260       if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1261
1262       if(debugStreaming){
1263         for (iTime = 0; iTime < fTimeTotal; iTime++)
1264           (*fDebugStream) << "TailCancellation"
1265                           << "col="  << iCol
1266                           << "row="  << iRow
1267                           << "\n";
1268       }
1269
1270       // Apply the tail cancelation via the digital filter
1271       DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1272
1273     } // while irow icol
1274
1275   return;
1276
1277 }
1278
1279 //_____________________________________________________________________________
1280 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1281 {
1282   //
1283   // Tail cancellation by deconvolution for PASA v4 TRF
1284   //
1285
1286   Float_t rates[2];
1287   Float_t coefficients[2];
1288
1289   // Initialization (coefficient = alpha, rates = lambda)
1290   Float_t r1 = 1.0;
1291   Float_t r2 = 1.0;
1292   Float_t c1 = 0.5;
1293   Float_t c2 = 0.5;
1294
1295   if (nexp == 1) {   // 1 Exponentials
1296     r1 = 1.156;
1297     r2 = 0.130;
1298     c1 = 0.066;
1299     c2 = 0.000;
1300   }
1301   if (nexp == 2) {   // 2 Exponentials
1302     Double_t par[4];
1303     fReconstructor->GetRecoParam()->GetTCParams(par);
1304     r1 = par[0];//1.156;
1305     r2 = par[1];//0.130;
1306     c1 = par[2];//0.114;
1307     c2 = par[3];//0.624;
1308   }
1309
1310   coefficients[0] = c1;
1311   coefficients[1] = c2;
1312
1313   Double_t dt = 0.1;
1314
1315   rates[0] = TMath::Exp(-dt/(r1));
1316   rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
1317
1318   Float_t reminder[2] = { .0, .0 };
1319   Float_t correction = 0.0;
1320   Float_t result     = 0.0;
1321
1322   for (int i = 0; i < nTime; i++) {
1323
1324     result = arr[i] - correction - fBaseline;    // No rescaling
1325     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1326
1327     correction = 0.0;
1328     for (int k = 0; k < 2; k++) {
1329       correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1330     }
1331
1332   }
1333
1334 }
1335
1336 //_____________________________________________________________________________
1337 void AliTRDclusterizer::ResetRecPoints()
1338 {
1339   //
1340   // Resets the list of rec points
1341   //
1342
1343   if (fRecPoints) {
1344     fRecPoints->Clear();
1345     fNoOfClusters = 0;
1346     //    delete fRecPoints;
1347   }
1348 }
1349
1350 //_____________________________________________________________________________
1351 TClonesArray *AliTRDclusterizer::RecPoints()
1352 {
1353   //
1354   // Returns the list of rec points
1355   //
1356
1357   if (!fRecPoints) {
1358     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1359       // determine number of clusters which has to be allocated
1360       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1361
1362       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1363     }
1364     //SetClustersOwner(kTRUE);
1365     AliTRDReconstructor::SetClusters(0x0);
1366   }
1367   return fRecPoints;
1368
1369 }
1370
1371 //_____________________________________________________________________________
1372 TClonesArray *AliTRDclusterizer::TrackletsArray()
1373 {
1374   //
1375   // Returns the list of rec points
1376   //
1377
1378   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1379     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1380     //SetClustersOwner(kTRUE);
1381     //AliTRDReconstructor::SetTracklets(0x0);
1382   }
1383   return fTracklets;
1384
1385 }
1386