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