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