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         if(Signals[0]<0)Signals[0]=0;
 900         if(Signals[2]<0)Signals[2]=0;
 901         Float_t  noiseSumThresh = fMinLeftRightCutSigma
 902           * fCalNoiseDetValue
 903           * fCalNoiseROC->GetValue(Max.col, Max.row);
 904         if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
 905         padStatus = 0;
 906         return kTRUE;
 907       }
 908     }
 909   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
 910     if(Signals[0]<0)Signals[0]=0;
 911     if(Signals[2]<0)Signals[2]=0;
 912     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
 913       Signals[2]=0;
 914       SetPadStatus(status[2], padStatus);
 915       return kTRUE;
 916     }
 917     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
 918       Signals[0]=0;
 919       SetPadStatus(status[0], padStatus);
 920       return kTRUE;
 921     }
 922     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
 923       Signals[1] = (Short_t)(fMaxThresh + 0.5f);
 924       SetPadStatus(status[1], padStatus);
 925       return kTRUE;
 926     }
 927   }
 928   return kFALSE;
 929 }
 930
 931 //_____________________________________________________________________________
 932 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
 933 {
 934   //
 935   // Look for 5 pad cluster with minimum in the middle
 936   // Gives back the ratio
 937   //
 938
 939   if (ThisMax.col >= fColMax - 3) return kFALSE;
 940   Float_t gain;
 941   if (ThisMax.col < fColMax - 5){
 942     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
 943     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
 944       return kFALSE;
 945   }
 946   if (ThisMax.col > 1) {
 947     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
 948     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
 949       return kFALSE;
 950   }
 951
 952   const Float_t kEpsilon = 0.01;
 953   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
 954       NeighbourMax.signals[1], NeighbourMax.signals[2]};
 955
 956   // Unfold the two maxima and set the signal on
 957   // the overlapping pad to the ratio
 958   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
 959   ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
 960   NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
 961   ThisMax.fivePad=kTRUE;
 962   NeighbourMax.fivePad=kTRUE;
 963   return kTRUE;
 964
 965 }
 966
 967 //_____________________________________________________________________________
 968 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 969 {
 970   //
 971   // Creates a cluster at the given position and saves it in fRecPoints
 972   //
 973
 974   Int_t nPadCount = 1;
 975   Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
 976   if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
 977
 978   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
 979   cluster.SetNPads(nPadCount);
 980   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
 981   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
 982   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
 983
 984   cluster.SetFivePad(Max.fivePad);
 985   // set pads status for the cluster
 986   UChar_t maskPosition = GetCorruption(Max.padStatus);
 987   if (maskPosition) {
 988     cluster.SetPadMaskedPosition(maskPosition);
 989     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
 990   }
 991
 992   // Transform the local cluster coordinates into calibrated
 993   // space point positions defined in the local tracking system.
 994   // Here the calibration for T0, Vdrift and ExB is applied as well.
 995   if(!fTransform->Transform(&cluster)) return;
 996   // Temporarily store the Max.Row, column and time bin of the center pad
 997   // Used to later on assign the track indices
 998   cluster.SetLabel(Max.row, 0);
 999   cluster.SetLabel(Max.col, 1);
1000   cluster.SetLabel(Max.time,2);
1001
1002   //needed for HLT reconstruction
1003   AddClusterToArray(&cluster);
1004
1005   // Store the index of the first cluster in the current ROC
1006   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1007
1008   fNoOfClusters++;
1009   fClusterROC++;
1010 }
1011
1012 //_____________________________________________________________________________
1013 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1014 {
1015   // Look to the right
1016   Int_t ii = 1;
1017   while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
1018     nPadCount++;
1019     ii++;
1020     if (Max.col < ii) break;
1021   }
1022   // Look to the left
1023   ii = 1;
1024   while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
1025     nPadCount++;
1026     ii++;
1027     if (Max.col+ii >= fColMax) break;
1028   }
1029
1030   // Store the amplitudes of the pads in the cluster for later analysis
1031   // and check whether one of these pads is masked in the database
1032   signals[2]=Max.signals[0];
1033   signals[3]=Max.signals[1];
1034   signals[4]=Max.signals[2];
1035   Float_t gain;
1036   for(Int_t i = 0; i<2; i++)
1037     {
1038       if(Max.col+i >= 3){
1039         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1040         signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1041       }
1042       if(Max.col+3-i < fColMax){
1043         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1044         signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1045       }
1046     }
1047   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1048     if ((jPad >= 0) && (jPad < fColMax))
1049       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1050       }*/
1051 }
1052
1053 //_____________________________________________________________________________
1054 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1055 {
1056   //
1057   // Add a cluster to the array
1058   //
1059
1060   Int_t n = RecPoints()->GetEntriesFast();
1061   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1062   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1063 }
1064
1065 //_____________________________________________________________________________
1066 void AliTRDclusterizer::AddTrackletsToArray()
1067 {
1068   //
1069   // Add the online tracklets of this chamber to the array
1070   //
1071
1072   UInt_t* trackletword;
1073   for(Int_t side=0; side<2; side++)
1074     {
1075       Int_t trkl=0;
1076       trackletword=fTrackletContainer[side];
1077       while(trackletword[trkl]>0){
1078         Int_t n = TrackletsArray()->GetEntriesFast();
1079         AliTRDtrackletWord tmp(trackletword[trkl]);
1080         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1081         trkl++;
1082       }
1083     }
1084 }
1085
1086 //_____________________________________________________________________________
1087 Bool_t AliTRDclusterizer::AddLabels()
1088 {
1089   //
1090   // Add the track indices to the found clusters
1091   //
1092
1093   const Int_t   kNclus  = 3;
1094   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1095   const Int_t   kNtrack = kNdict * kNclus;
1096
1097   Int_t iClusterROC = 0;
1098
1099   Int_t row  = 0;
1100   Int_t col  = 0;
1101   Int_t time = 0;
1102   Int_t iPad = 0;
1103
1104   // Temporary array to collect the track indices
1105   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1106
1107   // Loop through the dictionary arrays one-by-one
1108   // to keep memory consumption low
1109   AliTRDarrayDictionary *tracksIn = 0;  //mod
1110   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1111
1112     // tracksIn should be expanded beforehand!
1113     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1114
1115     // Loop though the clusters found in this ROC
1116     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1117
1118       AliTRDcluster *cluster = (AliTRDcluster *)
1119                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1120       row  = cluster->GetLabel(0);
1121       col  = cluster->GetLabel(1);
1122       time = cluster->GetLabel(2);
1123
1124       for (iPad = 0; iPad < kNclus; iPad++) {
1125         Int_t iPadCol = col - 1 + iPad;
1126         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1127         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1128       }
1129
1130     }
1131
1132   }
1133
1134   // Copy the track indices into the cluster
1135   // Loop though the clusters found in this ROC
1136   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1137
1138     AliTRDcluster *cluster = (AliTRDcluster *)
1139       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1140     cluster->SetLabel(-9999,0);
1141     cluster->SetLabel(-9999,1);
1142     cluster->SetLabel(-9999,2);
1143
1144     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1145
1146   }
1147
1148   delete [] idxTracks;
1149
1150   return kTRUE;
1151
1152 }
1153
1154 //_____________________________________________________________________________
1155 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1156 {
1157   //
1158   // Method to unfold neighbouring maxima.
1159   // The charge ratio on the overlapping pad is calculated
1160   // until there is no more change within the range given by eps.
1161   // The resulting ratio is then returned to the calling method.
1162   //
1163
1164   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1165   if (!calibration) {
1166     AliError("No AliTRDcalibDB instance available\n");
1167     return kFALSE;
1168   }
1169
1170   Int_t   irc                = 0;
1171   Int_t   itStep             = 0;                 // Count iteration steps
1172
1173   Double_t ratio             = 0.5;               // Start value for ratio
1174   Double_t prevRatio         = 0.0;               // Store previous ratio
1175
1176   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1177   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1178   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1179
1180   // Start the iteration
1181   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1182
1183     itStep++;
1184     prevRatio = ratio;
1185
1186     // Cluster position according to charge ratio
1187     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1188                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1189     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1190                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1191
1192     // Set cluster charge ratio
1193     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1194     Double_t ampLeft  = padSignal[1] / newSignal[1];
1195     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1196     Double_t ampRight = padSignal[3] / newSignal[1];
1197
1198     // Apply pad response to parameters
1199     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1200     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1201
1202     // Calculate new overlapping ratio
1203     ratio = TMath::Min((Double_t) 1.0
1204                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1205
1206   }
1207
1208   return ratio;
1209
1210 }
1211
1212 //_____________________________________________________________________________
1213 void AliTRDclusterizer::TailCancelation()
1214 {
1215   //
1216   // Applies the tail cancelation and gain factors:
1217   // Transform fDigits to fDigits
1218   //
1219
1220   Int_t iRow  = 0;
1221   Int_t iCol  = 0;
1222   Int_t iTime = 0;
1223
1224   Double_t *inADC = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1225   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1226
1227   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1228   Bool_t debugStreaming = fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1229   while(fIndexes->NextRCIndex(iRow, iCol))
1230     {
1231       Bool_t corrupted = kFALSE;
1232       if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1233
1234       // Save data into the temporary processing array and substract the baseline,
1235       // since DeConvExp does not expect a baseline
1236       for (iTime = 0; iTime < fTimeTotal; iTime++)
1237         inADC[iTime]   = fDigits->GetData(iRow,iCol,iTime)-fBaseline;
1238
1239       if(debugStreaming){
1240         for (iTime = 0; iTime < fTimeTotal; iTime++)
1241           (*fDebugStream) << "TailCancellation"
1242                           << "col="  << iCol
1243                           << "row="  << iRow
1244                           << "time=" << iTime
1245                           << "inADC=" << inADC[iTime]
1246                           << "outADC=" << outADC[iTime]
1247                           << "corrupted=" << corrupted
1248                           << "\n";
1249       }
1250
1251       if (!corrupted)
1252         {
1253           // Apply the tail cancelation via the digital filter
1254           // (only for non-coorupted pads)
1255           DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1256         }
1257       else memcpy(&outADC[0],&inADC[0],fTimeTotal*sizeof(inADC[0]));
1258
1259       // Save tailcancalled data and add the baseline
1260       for(iTime = 0; iTime < fTimeTotal; iTime++)
1261         fDigits->SetData(iRow,iCol,iTime,(Short_t)(outADC[iTime] + fBaseline + 0.5));
1262
1263     } // while irow icol
1264
1265   delete [] inADC;
1266   delete [] outADC;
1267
1268   return;
1269
1270 }
1271
1272 //_____________________________________________________________________________
1273 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1274                                   ,const Int_t n, const Int_t nexp)
1275 {
1276   //
1277   // Tail cancellation by deconvolution for PASA v4 TRF
1278   //
1279
1280   Double_t rates[2];
1281   Double_t coefficients[2];
1282
1283   // Initialization (coefficient = alpha, rates = lambda)
1284   Double_t r1 = 1.0;
1285   Double_t r2 = 1.0;
1286   Double_t c1 = 0.5;
1287   Double_t c2 = 0.5;
1288
1289   if (nexp == 1) {   // 1 Exponentials
1290     r1 = 1.156;
1291     r2 = 0.130;
1292     c1 = 0.066;
1293     c2 = 0.000;
1294   }
1295   if (nexp == 2) {   // 2 Exponentials
1296     Double_t par[4];
1297     fReconstructor->GetRecoParam()->GetTCParams(par);
1298     r1 = par[0];//1.156;
1299     r2 = par[1];//0.130;
1300     c1 = par[2];//0.114;
1301     c2 = par[3];//0.624;
1302   }
1303
1304   coefficients[0] = c1;
1305   coefficients[1] = c2;
1306
1307   Double_t dt = 0.1;
1308
1309   rates[0] = TMath::Exp(-dt/(r1));
1310   rates[1] = TMath::Exp(-dt/(r2));
1311
1312   Int_t i = 0;
1313   Int_t k = 0;
1314
1315   Double_t reminder[2];
1316   Double_t correction = 0.0;
1317   Double_t result     = 0.0;
1318
1319   // Attention: computation order is important
1320   for (k = 0; k < nexp; k++) {
1321     reminder[k] = 0.0;
1322   }
1323
1324   for (i = 0; i < n; i++) {
1325
1326     result    = (source[i] - correction);    // No rescaling
1327     target[i] = result;
1328
1329     for (k = 0; k < nexp; k++) {
1330       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1331     }
1332
1333     correction = 0.0;
1334     for (k = 0; k < nexp; k++) {
1335       correction += reminder[k];
1336     }
1337
1338   }
1339
1340 }
1341
1342 //_____________________________________________________________________________
1343 void AliTRDclusterizer::ResetRecPoints()
1344 {
1345   //
1346   // Resets the list of rec points
1347   //
1348
1349   if (fRecPoints) {
1350     fRecPoints->Delete();
1351     delete fRecPoints;
1352   }
1353 }
1354
1355 //_____________________________________________________________________________
1356 TClonesArray *AliTRDclusterizer::RecPoints()
1357 {
1358   //
1359   // Returns the list of rec points
1360   //
1361
1362   if (!fRecPoints) {
1363     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1364       // determine number of clusters which has to be allocated
1365       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1366
1367       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1368     }
1369     //SetClustersOwner(kTRUE);
1370     AliTRDReconstructor::SetClusters(0x0);
1371   }
1372   return fRecPoints;
1373
1374 }
1375
1376 //_____________________________________________________________________________
1377 TClonesArray *AliTRDclusterizer::TrackletsArray()
1378 {
1379   //
1380   // Returns the list of rec points
1381   //
1382
1383   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1384     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1385     //SetClustersOwner(kTRUE);
1386     //AliTRDReconstructor::SetTracklets(0x0);
1387   }
1388   return fTracklets;
1389
1390 }
1391