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