TRD/AliTRDclusterizer.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // TRD cluster finder                                                        //
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TF1.h>
  25 #include <TTree.h>
  26 #include <TH1.h>
  27 #include <TFile.h>
  28 #include <TClonesArray.h>
  29 #include <TObjArray.h>
  30
  31 #include "AliRunLoader.h"
  32 #include "AliLoader.h"
  33 #include "AliRawReader.h"
  34 #include "AliLog.h"
  35 #include "AliAlignObj.h"
  36
  37 #include "AliTRDclusterizer.h"
  38 #include "AliTRDcluster.h"
  39 #include "AliTRDReconstructor.h"
  40 #include "AliTRDgeometry.h"
  41 #include "AliTRDarraySignal.h"
  42 #include "AliTRDarrayDictionary.h"
  43 #include "AliTRDarrayADC.h"
  44 #include "AliTRDdigitsManager.h"
  45 #include "AliTRDrawData.h"
  46 #include "AliTRDcalibDB.h"
  47 #include "AliTRDrecoParam.h"
  48 #include "AliTRDCommonParam.h"
  49 #include "AliTRDtransform.h"
  50 #include "AliTRDSignalIndex.h"
  51 #include "AliTRDrawStreamBase.h"
  52 #include "AliTRDfeeParam.h"
  53
  54 #include "TTreeStream.h"
  55
  56 #include "Cal/AliTRDCalROC.h"
  57 #include "Cal/AliTRDCalDet.h"
  58 #include "Cal/AliTRDCalSingleChamberStatus.h"
  59
  60 ClassImp(AliTRDclusterizer)
  61
  62 //_____________________________________________________________________________
  63 AliTRDclusterizer::AliTRDclusterizer(AliTRDReconstructor *rec)
  64   :TNamed()
  65   ,fReconstructor(rec)
  66   ,fRunLoader(NULL)
  67   ,fClusterTree(NULL)
  68   ,fRecPoints(NULL)
  69   ,fTrackletTree(NULL)
  70   ,fDigitsManager(NULL)
  71   ,fTrackletContainer(NULL)
  72   ,fAddLabels(kTRUE)
  73   ,fRawVersion(2)
  74   ,fTransform(new AliTRDtransform(0))
  75   ,fLUTbin(0)
  76   ,fLUT(NULL)
  77   ,fDigitsIn(NULL)
  78   ,fIndexes(NULL)
  79   ,fADCthresh(0)
  80   ,fMaxThresh(0)
  81   ,fSigThresh(0)
  82   ,fMinMaxCutSigma(0)
  83   ,fMinLeftRightCutSigma(0)
  84   ,fLayer(0)
  85   ,fDet(0)
  86   ,fVolid(0)
  87   ,fColMax(0)
  88   ,fTimeTotal(0)
  89   ,fCalGainFactorROC(NULL)
  90   ,fCalGainFactorDetValue(0)
  91   ,fCalNoiseROC(NULL)
  92   ,fCalNoiseDetValue(0)
  93   ,fDigitsOut(NULL)
  94   ,fClusterROC(0)
  95   ,firstClusterROC(0)
  96 {
  97   //
  98   // AliTRDclusterizer default constructor
  99   //
 100
 101   AliTRDcalibDB *trd = 0x0;
 102   if (!(trd = AliTRDcalibDB::Instance())) {
 103     AliFatal("Could not get calibration object");
 104   }
 105
 106   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 107
 108   // Initialize debug stream
 109   if(fReconstructor){
 110     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 1){
 111       TDirectory *savedir = gDirectory;
 112       //fgGetDebugStream    = new TTreeSRedirector("TRD.ClusterizerDebug.root");
 113       savedir->cd();
 114     }
 115   }
 116
 117 }
 118
 119 //_____________________________________________________________________________
 120 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, AliTRDReconstructor *rec)
 121   :TNamed(name,title)
 122   ,fReconstructor(rec)
 123   ,fRunLoader(NULL)
 124   ,fClusterTree(NULL)
 125   ,fRecPoints(NULL)
 126   ,fTrackletTree(NULL)
 127   ,fDigitsManager(new AliTRDdigitsManager())
 128   ,fTrackletContainer(NULL)
 129   ,fAddLabels(kTRUE)
 130   ,fRawVersion(2)
 131   ,fTransform(new AliTRDtransform(0))
 132   ,fLUTbin(0)
 133   ,fLUT(NULL)
 134   ,fDigitsIn(NULL)
 135   ,fIndexes(NULL)
 136   ,fADCthresh(0)
 137   ,fMaxThresh(0)
 138   ,fSigThresh(0)
 139   ,fMinMaxCutSigma(0)
 140   ,fMinLeftRightCutSigma(0)
 141   ,fLayer(0)
 142   ,fDet(0)
 143   ,fVolid(0)
 144   ,fColMax(0)
 145   ,fTimeTotal(0)
 146   ,fCalGainFactorROC(NULL)
 147   ,fCalGainFactorDetValue(0)
 148   ,fCalNoiseROC(NULL)
 149   ,fCalNoiseDetValue(0)
 150   ,fDigitsOut(NULL)
 151   ,fClusterROC(0)
 152   ,firstClusterROC(0)
 153 {
 154   //
 155   // AliTRDclusterizer constructor
 156   //
 157
 158   AliTRDcalibDB *trd = 0x0;
 159   if (!(trd = AliTRDcalibDB::Instance())) {
 160     AliFatal("Could not get calibration object");
 161   }
 162
 163   fDigitsManager->CreateArrays();
 164
 165   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
 166
 167   FillLUT();
 168
 169 }
 170
 171 //_____________________________________________________________________________
 172 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
 173   :TNamed(c)
 174   ,fReconstructor(c.fReconstructor)
 175   ,fRunLoader(NULL)
 176   ,fClusterTree(NULL)
 177   ,fRecPoints(NULL)
 178   ,fTrackletTree(NULL)
 179   ,fDigitsManager(NULL)
 180   ,fTrackletContainer(NULL)
 181   ,fAddLabels(kTRUE)
 182   ,fRawVersion(2)
 183   ,fTransform(NULL)
 184   ,fLUTbin(0)
 185   ,fLUT(0)
 186   ,fDigitsIn(NULL)
 187   ,fIndexes(NULL)
 188   ,fADCthresh(0)
 189   ,fMaxThresh(0)
 190   ,fSigThresh(0)
 191   ,fMinMaxCutSigma(0)
 192   ,fMinLeftRightCutSigma(0)
 193   ,fLayer(0)
 194   ,fDet(0)
 195   ,fVolid(0)
 196   ,fColMax(0)
 197   ,fTimeTotal(0)
 198   ,fCalGainFactorROC(NULL)
 199   ,fCalGainFactorDetValue(0)
 200   ,fCalNoiseROC(NULL)
 201   ,fCalNoiseDetValue(0)
 202   ,fDigitsOut(NULL)
 203   ,fClusterROC(0)
 204   ,firstClusterROC(0)
 205 {
 206   //
 207   // AliTRDclusterizer copy constructor
 208   //
 209
 210   FillLUT();
 211
 212 }
 213
 214 //_____________________________________________________________________________
 215 AliTRDclusterizer::~AliTRDclusterizer()
 216 {
 217   //
 218   // AliTRDclusterizer destructor
 219   //
 220
 221   if (fRecPoints/* && IsClustersOwner()*/){
 222     fRecPoints->Delete();
 223     delete fRecPoints;
 224   }
 225
 226   if (fDigitsManager) {
 227     delete fDigitsManager;
 228     fDigitsManager = NULL;
 229   }
 230
 231   if (fTrackletContainer){
 232     delete fTrackletContainer;
 233     fTrackletContainer = NULL;
 234   }
 235
 236   if (fTransform){
 237     delete fTransform;
 238     fTransform     = NULL;
 239   }
 240
 241   if (fLUT) {
 242     delete [] fLUT;
 243     fLUT           = NULL;
 244   }
 245
 246   if (fDigitsOut) {
 247     delete fDigitsOut;
 248     fDigitsOut     = NULL;
 249   }
 250
 251 }
 252
 253 //_____________________________________________________________________________
 254 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
 255 {
 256   //
 257   // Assignment operator
 258   //
 259
 260   if (this != &c)
 261     {
 262       ((AliTRDclusterizer &) c).Copy(*this);
 263     }
 264
 265   return *this;
 266
 267 }
 268
 269 //_____________________________________________________________________________
 270 void AliTRDclusterizer::Copy(TObject &c) const
 271 {
 272   //
 273   // Copy function
 274   //
 275
 276   ((AliTRDclusterizer &) c).fClusterTree   = NULL;
 277   ((AliTRDclusterizer &) c).fRecPoints     = NULL;
 278   ((AliTRDclusterizer &) c).fTrackletTree  = NULL;
 279   ((AliTRDclusterizer &) c).fDigitsManager = NULL;
 280   ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
 281   ((AliTRDclusterizer &) c).fAddLabels     = fAddLabels;
 282   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
 283   ((AliTRDclusterizer &) c).fTransform     = NULL;
 284   ((AliTRDclusterizer &) c).fLUTbin        = 0;
 285   ((AliTRDclusterizer &) c).fLUT           = NULL;
 286   ((AliTRDclusterizer &) c).fDigitsIn      = NULL;
 287   ((AliTRDclusterizer &) c).fIndexes       = NULL;
 288   ((AliTRDclusterizer &) c).fADCthresh     = 0;
 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).fDigitsOut     = NULL;
 303   ((AliTRDclusterizer &) c).fClusterROC    = 0;
 304   ((AliTRDclusterizer &) c).firstClusterROC= 0;
 305
 306 }
 307
 308 //_____________________________________________________________________________
 309 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
 310 {
 311   //
 312   // Opens the AliROOT file. Output and input are in the same file
 313   //
 314
 315   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 316   fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 317
 318   if (!fRunLoader) {
 319     fRunLoader = AliRunLoader::Open(name);
 320   }
 321
 322   if (!fRunLoader) {
 323     AliError(Form("Can not open session for file %s.",name));
 324     return kFALSE;
 325   }
 326
 327   OpenInput(nEvent);
 328   OpenOutput();
 329
 330   return kTRUE;
 331
 332 }
 333
 334 //_____________________________________________________________________________
 335 Bool_t AliTRDclusterizer::OpenOutput()
 336 {
 337   //
 338   // Open the output file
 339   //
 340
 341   if (!fReconstructor->IsWritingClusters()) return kTRUE;
 342
 343   TObjArray *ioArray = 0x0;
 344
 345   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 346   loader->MakeTree("R");
 347
 348   fClusterTree = loader->TreeR();
 349   fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 350
 351   return kTRUE;
 352
 353 }
 354
 355 //_____________________________________________________________________________
 356 Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
 357 {
 358   //
 359   // Connect the output tree
 360   //
 361
 362   // clusters writing
 363   if (fReconstructor->IsWritingClusters()){
 364     TObjArray *ioArray = 0x0;
 365     fClusterTree = clusterTree;
 366     fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
 367   }
 368
 369   // tracklet writing
 370   if (fReconstructor->IsWritingTracklets()){
 371     TString evfoldname = AliConfig::GetDefaultEventFolderName();
 372     fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
 373
 374     if (!fRunLoader) {
 375       fRunLoader = AliRunLoader::Open("galice.root");
 376     }
 377     if (!fRunLoader) {
 378       AliError(Form("Can not open session for file galice.root."));
 379       return kFALSE;
 380     }
 381
 382     UInt_t **leaves = new UInt_t *[2];
 383     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 384     if (!dl) {
 385       AliError("Could not get the tracklets data loader!");
 386       dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
 387       fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
 388     }
 389     else {
 390       fTrackletTree = dl->Tree();
 391       if (!fTrackletTree)
 392         {
 393         dl->MakeTree();
 394         fTrackletTree = dl->Tree();
 395         }
 396       TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 397       if (!trkbranch)
 398         fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 399     }
 400   }
 401
 402   return kTRUE;
 403
 404 }
 405
 406 //_____________________________________________________________________________
 407 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
 408 {
 409   //
 410   // Opens a ROOT-file with TRD-hits and reads in the digits-tree
 411   //
 412
 413   // Import the Trees for the event nEvent in the file
 414   fRunLoader->GetEvent(nEvent);
 415
 416   return kTRUE;
 417
 418 }
 419
 420 //_____________________________________________________________________________
 421 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
 422 {
 423   //
 424   // Fills TRDcluster branch in the tree with the clusters
 425   // found in detector = det. For det=-1 writes the tree.
 426   //
 427
 428   if ((det <                      -1) ||
 429       (det >= AliTRDgeometry::Ndet())) {
 430     AliError(Form("Unexpected detector index %d.\n",det));
 431     return kFALSE;
 432   }
 433
 434   TObjArray *ioArray = new TObjArray(400);
 435   TBranch *branch = fClusterTree->GetBranch("TRDcluster");
 436   if (!branch) {
 437     branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
 438   } else branch->SetAddress(&ioArray);
 439
 440   Int_t nRecPoints = RecPoints()->GetEntriesFast();
 441   if(det >= 0){
 442     for (Int_t i = 0; i < nRecPoints; i++) {
 443       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 444       if(det != c->GetDetector()) continue;
 445       ioArray->AddLast(c);
 446     }
 447     fClusterTree->Fill();
 448   } else {
 449
 450     Int_t detOld = -1;
 451     for (Int_t i = 0; i < nRecPoints; i++) {
 452       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
 453       if(c->GetDetector() != detOld){
 454         fClusterTree->Fill();
 455         ioArray->Clear();
 456         detOld = c->GetDetector();
 457       }
 458       ioArray->AddLast(c);
 459     }
 460   }
 461   delete ioArray;
 462
 463   return kTRUE;
 464
 465 }
 466
 467 //_____________________________________________________________________________
 468 Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
 469 {
 470   //
 471   // Write the raw data tracklets into seperate file
 472   //
 473
 474   UInt_t **leaves = new UInt_t *[2];
 475   for (Int_t i=0; i<2 ;i++){
 476     leaves[i] = new UInt_t[258];
 477     leaves[i][0] = det; // det
 478     leaves[i][1] = i;   // side
 479     memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
 480   }
 481
 482   if (!fTrackletTree){
 483     AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 484     dl->MakeTree();
 485     fTrackletTree = dl->Tree();
 486   }
 487
 488   TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
 489   if (!trkbranch) {
 490     trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
 491   }
 492
 493   for (Int_t i=0; i<2; i++){
 494     if (leaves[i][2]>0) {
 495       trkbranch->SetAddress(leaves[i]);
 496       fTrackletTree->Fill();
 497     }
 498   }
 499
 500   AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
 501   dl->WriteData("OVERWRITE");
 502   //dl->Unload();
 503   delete [] leaves;
 504
 505   return kTRUE;
 506
 507 }
 508
 509 //_____________________________________________________________________________
 510 Bool_t AliTRDclusterizer::ReadDigits()
 511 {
 512   //
 513   // Reads the digits arrays from the input aliroot file
 514   //
 515
 516   if (!fRunLoader) {
 517     AliError("No run loader available");
 518     return kFALSE;
 519   }
 520
 521   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 522   if (!loader->TreeD()) {
 523     loader->LoadDigits();
 524   }
 525
 526   // Read in the digit arrays
 527   return (fDigitsManager->ReadDigits(loader->TreeD()));
 528
 529 }
 530
 531 //_____________________________________________________________________________
 532 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
 533 {
 534   //
 535   // Reads the digits arrays from the input tree
 536   //
 537
 538   // Read in the digit arrays
 539   return (fDigitsManager->ReadDigits(digitsTree));
 540
 541 }
 542
 543 //_____________________________________________________________________________
 544 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
 545 {
 546   //
 547   // Reads the digits arrays from the ddl file
 548   //
 549
 550   AliTRDrawData raw;
 551   fDigitsManager = raw.Raw2Digits(rawReader);
 552
 553   return kTRUE;
 554
 555 }
 556
 557 //_____________________________________________________________________________
 558 Bool_t AliTRDclusterizer::MakeClusters()
 559 {
 560   //
 561   // Creates clusters from digits
 562   //
 563
 564   // Propagate info from the digits manager
 565   if (fAddLabels == kTRUE){
 566     fAddLabels = fDigitsManager->UsesDictionaries();
 567   }
 568
 569   Bool_t fReturn = kTRUE;
 570   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
 571
 572     AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
 573     // This is to take care of switched off super modules
 574     if (!digitsIn->HasData()) continue;
 575     digitsIn->Expand();
 576     digitsIn->DeleteNegatives();  // Restore digits array to >=0 values
 577     AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 578     if (indexes->IsAllocated() == kFALSE){
 579       fDigitsManager->BuildIndexes(i);
 580     }
 581
 582     Bool_t fR = kFALSE;
 583     if (indexes->HasEntry()){
 584       if (fAddLabels){
 585         for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
 586           AliTRDarrayDictionary *tracksIn = 0; //mod
 587           tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict);  //mod
 588           tracksIn->Expand();
 589         }
 590       }
 591       fR = MakeClusters(i);
 592       fReturn = fR && fReturn;
 593     }
 594
 595     //if (fR == kFALSE){
 596     //  if(IsWritingClusters()) WriteClusters(i);
 597     //  ResetRecPoints();
 598     //}
 599
 600     // No compress just remove
 601     fDigitsManager->RemoveDigits(i);
 602     fDigitsManager->RemoveDictionaries(i);
 603     fDigitsManager->ClearIndexes(i);
 604   }
 605
 606   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 607
 608   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
 609
 610   return fReturn;
 611
 612 }
 613
 614 //_____________________________________________________________________________
 615 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
 616 {
 617   //
 618   // Creates clusters from raw data
 619   //
 620
 621   return Raw2ClustersChamber(rawReader);
 622
 623 }
 624
 625 //_____________________________________________________________________________
 626 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 627 {
 628   //
 629   // Creates clusters from raw data
 630   //
 631
 632   // Create the digits manager
 633   if (!fDigitsManager){
 634     fDigitsManager = new AliTRDdigitsManager();
 635     fDigitsManager->CreateArrays();
 636   }
 637
 638   fDigitsManager->SetUseDictionaries(fAddLabels);
 639
 640   // tracklet container for raw tracklet writing
 641   if (!fTrackletContainer && fReconstructor->IsWritingTracklets()) {
 642     // maximum tracklets for one HC
 643     const Int_t kTrackletChmb=256;
 644     fTrackletContainer = new UInt_t *[2];
 645     fTrackletContainer[0] = new UInt_t[kTrackletChmb];
 646     fTrackletContainer[1] = new UInt_t[kTrackletChmb];
 647   }
 648
 649
 650   AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader);
 651
 652   AliInfo(Form("Stream version: %s", input->IsA()->GetName()));
 653
 654   Int_t det    = 0;
 655   while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
 656     Bool_t iclusterBranch = kFALSE;
 657     if (fDigitsManager->GetIndexes(det)->HasEntry()){
 658     iclusterBranch = MakeClusters(det);
 659     }
 660
 661     fDigitsManager->RemoveDigits(det);
 662     fDigitsManager->RemoveDictionaries(det);
 663     fDigitsManager->ClearIndexes(det);
 664
 665     if (!fReconstructor->IsWritingTracklets()) continue;
 666     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
 667   }
 668   if (fReconstructor->IsWritingTracklets()){
 669     delete [] fTrackletContainer[0];
 670     delete [] fTrackletContainer[1];
 671     delete [] fTrackletContainer;
 672     fTrackletContainer = NULL;
 673   }
 674
 675   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 676
 677   delete fDigitsManager;
 678   fDigitsManager = NULL;
 679
 680   delete input;
 681   input = NULL;
 682
 683   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
 684   return kTRUE;
 685
 686 }
 687
 688 //_____________________________________________________________________________
 689 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
 690 {
 691   //
 692   // Check if a pad is masked
 693   //
 694
 695   UChar_t status = 0;
 696
 697   if(signal>0 && TESTBIT(signal, 10)){
 698     CLRBIT(signal, 10);
 699     for(int ibit=0; ibit<4; ibit++){
 700       if(TESTBIT(signal, 11+ibit)){
 701         SETBIT(status, ibit);
 702         CLRBIT(signal, 11+ibit);
 703       }
 704     }
 705   }
 706   return status;
 707 }
 708
 709 //_____________________________________________________________________________
 710 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out){
 711   //
 712   // Set the pad status into out
 713   // First three bits are needed for the position encoding
 714   //
 715   out |= status << 3;
 716 }
 717
 718 //_____________________________________________________________________________
 719 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
 720   //
 721   // return the staus encoding of the corrupted pad
 722   //
 723   return static_cast<UChar_t>(encoding >> 3);
 724 }
 725
 726 //_____________________________________________________________________________
 727 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
 728   //
 729   // Return the position of the corruption
 730   //
 731   return encoding & 7;
 732 }
 733
 734 //_____________________________________________________________________________
 735 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 736 {
 737   //
 738   // Generates the cluster.
 739   //
 740
 741   // Get the digits
 742   //   digits should be expanded beforehand!
 743   //   digitsIn->Expand();
 744   fDigitsIn = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
 745
 746   // This is to take care of switched off super modules
 747   if (!fDigitsIn->HasData())
 748     {
 749       return kFALSE;
 750     }
 751
 752   fIndexes = fDigitsManager->GetIndexes(det);
 753   if (fIndexes->IsAllocated() == kFALSE)
 754     {
 755       AliError("Indexes do not exist!");
 756       return kFALSE;
 757     }
 758
 759   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
 760   if (!calibration)
 761     {
 762       AliFatal("No AliTRDcalibDB instance available\n");
 763       return kFALSE;
 764     }
 765
 766   fADCthresh = 0;
 767
 768   if (!fReconstructor){
 769     AliError("Reconstructor not set\n");
 770     return kFALSE;
 771   }
 772
 773   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDReconstructor::kClusterizer);
 774
 775   fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
 776   fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
 777   fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
 778   fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
 779
 780   Int_t istack  = fIndexes->GetStack();
 781   fLayer  = fIndexes->GetLayer();
 782   Int_t isector = fIndexes->GetSM();
 783
 784   // Start clustering in the chamber
 785
 786   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
 787   if (fDet != det) {
 788     AliError("Strange Detector number Missmatch!");
 789     return kFALSE;
 790   }
 791
 792   // TRD space point transformation
 793   fTransform->SetDetector(det);
 794
 795   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
 796   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
 797   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
 798
 799   fColMax    = fDigitsIn->GetNcol();
 800   Int_t nRowMax    = fDigitsIn->GetNrow();
 801   fTimeTotal = fDigitsIn->GetNtime();
 802
 803   // Detector wise calibration object for the gain factors
 804   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
 805   // Calibration object with pad wise values for the gain factors
 806   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
 807   // Calibration value for chamber wise gain factor
 808   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
 809
 810   // Detector wise calibration object for the noise
 811   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
 812   // Calibration object with pad wise values for the noise
 813   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
 814   // Calibration value for chamber wise noise
 815   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
 816
 817   if(fDigitsOut) delete fDigitsOut;
 818   fDigitsOut = new AliTRDarraySignal(nRowMax, fColMax, fTimeTotal);
 819
 820   firstClusterROC = -1;
 821   fClusterROC     =  0;
 822
 823   // Apply the gain and the tail cancelation via digital filter
 824   TailCancelation();
 825
 826   ClusterizerStruct curr, last;
 827   last.Row = -1;
 828   Int_t nMaximas = 0, nCorrupted = 0;
 829   Double_t Ratio = 1;
 830
 831   // Here the clusterfining is happening
 832
 833   for(curr.Time = 0; curr.Time < fTimeTotal; curr.Time++)
 834     while(fIndexes->NextRCIndex(curr.Row, curr.Col))
 835       if(IsMaximum(curr.Row, curr.Col, curr.Time, curr.padStatus, &curr.Signals[0]))
 836         {
 837           if(last.Row>-1)
 838                 {
 839                   last.Signals[0] *= Ratio;
 840                   if(curr.Row==last.Row && curr.Col==last.Col+2)
 841                     {
 842                       if(IsFivePadCluster(last.Row, last.Col, last.Time, &last.Signals[0], &curr.Signals[0], Ratio))
 843                         {
 844                           last.Signals[2] *= Ratio;
 845                           Ratio = 1 - Ratio;
 846                         }else Ratio = 1;
 847                     }else Ratio = 1;
 848                   CreateCluster(last.Row, last.Col, last.Time, &last.Signals[0], last.padStatus);
 849                 }
 850           last=curr;
 851         }
 852   if(last.Row>-1)
 853     {
 854       last.Signals[0] *= Ratio;
 855       CreateCluster(last.Row, last.Col, last.Time, &last.Signals[0], last.padStatus);
 856     }
 857
 858   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 2){
 859     (*fDebugStream) << "MakeClusters"
 860                     << "Detector="   << det
 861                     << "NMaxima="    << nMaximas
 862                     << "NClusters="  << fClusterROC
 863                     << "NCorrupted=" << nCorrupted
 864                     << "\n";
 865   }
 866
 867   if (fAddLabels) {
 868     AddLabels(fDet,firstClusterROC,fClusterROC);
 869   }
 870
 871   return kTRUE;
 872
 873 }
 874
 875 //_____________________________________________________________________________
 876 Bool_t AliTRDclusterizer::IsMaximum(const Int_t row, const Int_t col, const Int_t time,
 877                                     UChar_t &pasStatus, Double_t *const Signals)
 878 {
 879   //
 880   // Returns true if this row,col,time combination is a maximum.
 881   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
 882   //
 883
 884   Signals[1] = fDigitsOut->GetData(row,col,time);
 885   if(Signals[1] < fMaxThresh) return kFALSE;
 886
 887   Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(col,row);
 888   if (Signals[1] < noiseMiddleThresh) return kFALSE;
 889
 890   if (col + 1 >= fColMax || col < 1) return kFALSE;
 891   UChar_t status[3]={0, 0, 0};
 892   pasStatus = 0;
 893
 894   status[1] = fDigitsIn->GetPadStatus(row,col,time);
 895   //if(status[1]) SETBIT(pasStatus, AliTRDcluster::kMaskedCenter);//TR: mod: this is already done by SetPadStatus
 896
 897   Signals[2] = fDigitsOut->GetData(row,col+1,time);
 898   status[2] = fDigitsIn->GetPadStatus(row,col+1,time);
 899   //if(status[2]) SETBIT(pasStatus, AliTRDcluster::kMaskedLeft);//TR: mod: this is already done by SetPadStatus
 900
 901   Signals[0] = fDigitsOut->GetData(row,col-1,time);
 902   status[0] = fDigitsIn->GetPadStatus(row,col-1,time);
 903   //if(status[0]) SETBIT(pasStatus, AliTRDcluster::kMaskedRight);//TR: mod: this is already done by SetPadStatus
 904
 905   // reject candidates with more than 1 problematic pad
 906   if(pasStatus >= 3) return kFALSE;
 907
 908   if (!status[1]) { // good central pad
 909     if (!pasStatus) { // all pads are OK
 910       if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
 911         if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
 912           Float_t  noiseSumThresh = fMinLeftRightCutSigma
 913             * fCalNoiseDetValue
 914             * fCalNoiseROC->GetValue(col,row);
 915           if ((Signals[2]+Signals[0]+Signals[1]) >= noiseSumThresh)
 916             return kTRUE;
 917         }
 918       }
 919     } else { // one of the neighbouring pads are bad
 920       if (status[2] && Signals[0] < Signals[1] && Signals[0] >= fSigThresh) {
 921         fDigitsOut->SetData(row, col+1, time, 0.);//TR: mod: was: SetData(row, col, time+1, 0.)
 922         SetPadStatus(status[2], pasStatus);
 923         return kTRUE;
 924       }
 925       else if (status[0] && Signals[2] <= Signals[1] && Signals[2] >= fSigThresh) {
 926         fDigitsOut->SetData(row, col-1, time, 0.);//TR: mod: was: SetData(row, col, time-1, 0.)
 927         SetPadStatus(status[0], pasStatus);
 928         return kTRUE;
 929       }
 930     }
 931   }
 932   else { // wrong maximum pad
 933     if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
 934       fDigitsOut->SetData(row,col,time,fMaxThresh);
 935       SetPadStatus(status[1], pasStatus);
 936       return kTRUE;
 937     }
 938   }
 939   return kFALSE;
 940 }
 941
 942 //_____________________________________________________________________________
 943 Bool_t AliTRDclusterizer::IsFivePadCluster(const Int_t row, const Int_t col, const Int_t time,
 944                                            Double_t *SignalsThisMax, Double_t *SignalsNeighbourMax, Double_t &ratio)
 945 {
 946   //
 947   // Look for 5 pad cluster with minimum in the middle
 948   // Gives back the ratio
 949   //
 950
 951   if (col < fColMax - 3){
 952     if (col < fColMax - 5){
 953       if (fDigitsOut->GetData(row,col+4,time) >= fSigThresh)
 954         return kFALSE;
 955     }
 956     if (col > 1) {
 957       if (fDigitsOut->GetData(row,col-2,time) >= fSigThresh)
 958         return kFALSE;
 959         }
 960
 961   //if (fSignalsThisMax[1] >= 0){ //TR: mod
 962
 963     const Float_t kEpsilon = 0.01;
 964     Double_t padSignal[5] = {SignalsThisMax[0],SignalsThisMax[1],SignalsThisMax[2],
 965                              SignalsNeighbourMax[1], SignalsNeighbourMax[2]};
 966
 967     // Unfold the two maxima and set the signal on
 968     // the overlapping pad to the ratio
 969     ratio = Unfold(kEpsilon,fLayer,padSignal);
 970     return kTRUE;
 971   }
 972   return kFALSE;
 973 }
 974
 975 //_____________________________________________________________________________
 976 void AliTRDclusterizer::CreateCluster(const Int_t row, const Int_t col, const Int_t time,
 977                                         const Double_t* const clusterSignal, const UChar_t pasStatus)
 978 {
 979   //
 980   // Creates a cluster at the given position and saves it in fRecPoint
 981   //
 982
 983   const Int_t kNsig = 5;
 984   Double_t padSignal[kNsig];
 985
 986   // The position of the cluster in COL direction relative to the center pad (pad units)
 987   Double_t clusterPosCol = 0.0;
 988   if (fReconstructor->GetRecoParam()->IsLUT()) {
 989     // Calculate the position of the cluster by using the
 990     // lookup table method
 991     clusterPosCol = LUTposition(fLayer,clusterSignal[0]
 992                                 ,clusterSignal[1]
 993                                 ,clusterSignal[2]);
 994   }
 995   else {
 996     // Calculate the position of the cluster by using the
 997     // center of gravity method
 998     padSignal[1] = clusterSignal[0];
 999     padSignal[2] = clusterSignal[1];
1000     padSignal[3] = clusterSignal[2];
1001     if(col > 2){
1002       padSignal[0] = fDigitsOut->GetData(row,col-2,time);
1003       if(padSignal[0]>= padSignal[1])
1004         padSignal[0] = 0;
1005     }
1006     if(col < fColMax - 3){
1007       padSignal[4] = fDigitsOut->GetData(row,col+2,time);
1008       if(padSignal[4]>= padSignal[3])
1009         padSignal[4] = 0;
1010     }
1011     clusterPosCol = GetCOG(padSignal);
1012   }
1013
1014   // Count the number of pads in the cluster
1015   Int_t nPadCount = 1;
1016   // Look to the right
1017   Int_t ii = 1;
1018   while (fDigitsOut->GetData(row, col-ii, time) >= fSigThresh) {
1019     nPadCount++;
1020     ii++;
1021     if (col-ii < 0) break;
1022   }
1023   // Look to the left
1024   ii = 1;
1025   while (fDigitsOut->GetData(row, col+ii, time) >= fSigThresh) {
1026     nPadCount++;
1027     ii++;
1028     if (col+ii >= fColMax) break;
1029   }
1030
1031   // Store the amplitudes of the pads in the cluster for later analysis
1032   // and check whether one of these pads is masked in the database
1033   Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
1034   for(Int_t i = 0; i++; i<3)
1035     signals[i+2] = TMath::Nint(clusterSignal[i]);
1036   for(Int_t i = 0; i++; i<2)
1037     {
1038       if(col+i >= 3)
1039         signals[i] = TMath::Nint(fDigitsOut->GetData(row,col-3+i,time));
1040       if(col+3-i < fColMax)
1041         signals[7-i] = TMath::Nint(fDigitsOut->GetData(row,col+3-i,time));
1042     }
1043   /*for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
1044     if ((jPad >= 0) && (jPad < fColMax))
1045       signals[jPad-col+3] = TMath::Nint(fDigitsOut->GetData(row,jPad,time));
1046       }*/
1047
1048   // Transform the local cluster coordinates into calibrated
1049   // space point positions defined in the local tracking system.
1050   // Here the calibration for T0, Vdrift and ExB is applied as well.
1051   Double_t clusterXYZ[6];
1052   clusterXYZ[0] = clusterPosCol;
1053   clusterXYZ[1] = clusterSignal[2];
1054   clusterXYZ[2] = clusterSignal[1];
1055   clusterXYZ[3] = clusterSignal[0];
1056   clusterXYZ[4] = 0.0;
1057   clusterXYZ[5] = 0.0;
1058   Int_t    clusterRCT[3];
1059   clusterRCT[0] = row;
1060   clusterRCT[1] = col;
1061   clusterRCT[2] = 0;
1062
1063   Bool_t out = kTRUE;
1064   if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),out,0)) {
1065
1066     // Add the cluster to the output array
1067     // The track indices will be stored later
1068     Float_t clusterPos[3];
1069     clusterPos[0] = clusterXYZ[0];
1070     clusterPos[1] = clusterXYZ[1];
1071     clusterPos[2] = clusterXYZ[2];
1072     Float_t clusterSig[2];
1073     clusterSig[0] = clusterXYZ[4];
1074     clusterSig[1] = clusterXYZ[5];
1075     Double_t clusterCharge  = clusterXYZ[3];
1076     Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
1077
1078     Int_t n = RecPoints()->GetEntriesFast();
1079     AliTRDcluster *cluster = new((*RecPoints())[n]) AliTRDcluster(
1080                                                                   fDet,
1081                                                                   clusterCharge, clusterPos, clusterSig,
1082                                                                   0x0,
1083                                                                   ((Char_t) nPadCount),
1084                                                                   signals,
1085                                                                   ((UChar_t) col), ((UChar_t) row), ((UChar_t) time),
1086                                                                   clusterTimeBin, clusterPosCol,
1087                                                                   fVolid);
1088     cluster->SetInChamber(!out);
1089
1090     UChar_t maskPosition = GetCorruption(pasStatus);
1091     UChar_t padstatus = GetPadStatus(pasStatus);
1092     if (maskPosition) {
1093       cluster->SetPadMaskedPosition(maskPosition);
1094       cluster->SetPadMaskedStatus(padstatus);
1095     }
1096
1097     // Temporarily store the row, column and time bin of the center pad
1098     // Used to later on assign the track indices
1099     cluster->SetLabel(row, 0);
1100     cluster->SetLabel(col, 1);
1101     cluster->SetLabel(time,2);
1102
1103     // Store the index of the first cluster in the current ROC
1104     if (firstClusterROC < 0) {
1105       firstClusterROC = RecPoints()->GetEntriesFast() - 1;
1106     }
1107
1108     // Count the number of cluster in the current ROC
1109     fClusterROC++;
1110   }
1111
1112 }
1113
1114 //_____________________________________________________________________________
1115 Bool_t AliTRDclusterizer::AddLabels(const Int_t idet, const Int_t firstClusterROC, const Int_t nClusterROC)
1116 {
1117   //
1118   // Add the track indices to the found clusters
1119   //
1120
1121   const Int_t   kNclus  = 3;
1122   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1123   const Int_t   kNtrack = kNdict * kNclus;
1124
1125   Int_t iClusterROC = 0;
1126
1127   Int_t row  = 0;
1128   Int_t col  = 0;
1129   Int_t time = 0;
1130   Int_t iPad = 0;
1131
1132   // Temporary array to collect the track indices
1133   Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1134
1135   // Loop through the dictionary arrays one-by-one
1136   // to keep memory consumption low
1137   AliTRDarrayDictionary *tracksIn = 0;  //mod
1138   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1139
1140     // tracksIn should be expanded beforehand!
1141     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(idet,iDict);
1142
1143     // Loop though the clusters found in this ROC
1144     for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1145
1146       AliTRDcluster *cluster = (AliTRDcluster *)
1147                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1148       row  = cluster->GetLabel(0);
1149       col  = cluster->GetLabel(1);
1150       time = cluster->GetLabel(2);
1151
1152       for (iPad = 0; iPad < kNclus; iPad++) {
1153         Int_t iPadCol = col - 1 + iPad;
1154         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1155         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1156       }
1157
1158     }
1159
1160   }
1161
1162   // Copy the track indices into the cluster
1163   // Loop though the clusters found in this ROC
1164   for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1165
1166     AliTRDcluster *cluster = (AliTRDcluster *)
1167       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1168     cluster->SetLabel(-9999,0);
1169     cluster->SetLabel(-9999,1);
1170     cluster->SetLabel(-9999,2);
1171
1172     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1173
1174   }
1175
1176   delete [] idxTracks;
1177
1178   return kTRUE;
1179
1180 }
1181
1182 //_____________________________________________________________________________
1183 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
1184 {
1185   //
1186   // Get COG position
1187   // Used for clusters with more than 3 pads - where LUT not applicable
1188   //
1189
1190   Double_t sum = signal[0]
1191               + signal[1]
1192               + signal[2]
1193               + signal[3]
1194               + signal[4];
1195
1196   // ???????????? CBL
1197   // Go to 3 pad COG ????
1198   // ???????????? CBL
1199   Double_t res = (0.0 * (-signal[0] + signal[4])
1200                       + (-signal[1] + signal[3])) / sum;
1201
1202   return res;
1203
1204 }
1205
1206 //_____________________________________________________________________________
1207 Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const
1208 {
1209   //
1210   // Method to unfold neighbouring maxima.
1211   // The charge ratio on the overlapping pad is calculated
1212   // until there is no more change within the range given by eps.
1213   // The resulting ratio is then returned to the calling method.
1214   //
1215
1216   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1217   if (!calibration) {
1218     AliError("No AliTRDcalibDB instance available\n");
1219     return kFALSE;
1220   }
1221
1222   Int_t   irc                = 0;
1223   Int_t   itStep             = 0;                 // Count iteration steps
1224
1225   Double_t ratio             = 0.5;               // Start value for ratio
1226   Double_t prevRatio         = 0.0;               // Store previous ratio
1227
1228   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1229   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1230   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1231
1232   // Start the iteration
1233   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1234
1235     itStep++;
1236     prevRatio = ratio;
1237
1238     // Cluster position according to charge ratio
1239     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1240                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1241     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1242                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1243
1244     // Set cluster charge ratio
1245     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1246     Double_t ampLeft  = padSignal[1] / newSignal[1];
1247     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1248     Double_t ampRight = padSignal[3] / newSignal[1];
1249
1250     // Apply pad response to parameters
1251     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1252     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1253
1254     // Calculate new overlapping ratio
1255     ratio = TMath::Min((Double_t) 1.0
1256                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1257
1258   }
1259
1260   return ratio;
1261
1262 }
1263
1264 //_____________________________________________________________________________
1265 void AliTRDclusterizer::TailCancelation()
1266 {
1267   //
1268   // Applies the tail cancelation and gain factors:
1269   // Transform fDigitsIn to fDigitsOut
1270   //
1271
1272   Int_t iRow  = 0;
1273   Int_t iCol  = 0;
1274   Int_t iTime = 0;
1275
1276   Double_t *inADC  = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1277   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1278   fIndexes->ResetCounters();
1279   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDReconstructor::kClusterizer);
1280   while(fIndexes->NextRCIndex(iRow, iCol))
1281     {
1282       Float_t  fCalGainFactorROCValue = fCalGainFactorROC->GetValue(iCol,iRow);
1283       Double_t gain                  = fCalGainFactorDetValue
1284                                     * fCalGainFactorROCValue;
1285
1286       Bool_t corrupted = kFALSE;
1287       for (iTime = 0; iTime < fTimeTotal; iTime++)
1288         {
1289           // Apply gain gain factor
1290           inADC[iTime]   = fDigitsIn->GetDataB(iRow,iCol,iTime);
1291           if (fDigitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
1292           inADC[iTime]  /= gain;
1293           outADC[iTime]  = inADC[iTime];
1294           if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 7){
1295             (*fDebugStream) << "TailCancellation"
1296                               << "col="  << iCol
1297                               << "row="  << iRow
1298                               << "time=" << iTime
1299                               << "inADC=" << inADC[iTime]
1300                               << "gain=" << gain
1301                               << "outADC=" << outADC[iTime]
1302                               << "corrupted=" << corrupted
1303                               << "\n";
1304           }
1305         }
1306       if (!corrupted)
1307         {
1308           // Apply the tail cancelation via the digital filter
1309           // (only for non-coorupted pads)
1310           if (fReconstructor->GetRecoParam()->IsTailCancelation())
1311             DeConvExp(inADC,outADC,fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1312         }
1313
1314       for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime))
1315         {
1316           // Store the amplitude of the digit if above threshold
1317           if (outADC[iTime] > fADCthresh)
1318             fDigitsOut->SetData(iRow,iCol,iTime,outADC[iTime]);
1319         } // while itime
1320
1321     } // while irow icol
1322
1323   delete [] inADC;
1324   delete [] outADC;
1325
1326   return;
1327
1328 }
1329
1330 //_____________________________________________________________________________
1331 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1332                                   ,const Int_t n, const Int_t nexp)
1333 {
1334   //
1335   // Tail cancellation by deconvolution for PASA v4 TRF
1336   //
1337
1338   Double_t rates[2];
1339   Double_t coefficients[2];
1340
1341   // Initialization (coefficient = alpha, rates = lambda)
1342   Double_t r1 = 1.0;
1343   Double_t r2 = 1.0;
1344   Double_t c1 = 0.5;
1345   Double_t c2 = 0.5;
1346
1347   if (nexp == 1) {   // 1 Exponentials
1348     r1 = 1.156;
1349     r2 = 0.130;
1350     c1 = 0.066;
1351     c2 = 0.000;
1352   }
1353   if (nexp == 2) {   // 2 Exponentials
1354     Double_t par[4];
1355     fReconstructor->GetTCParams(par);
1356     r1 = par[0];//1.156;
1357     r2 = par[1];//0.130;
1358     c1 = par[2];//0.114;
1359     c2 = par[3];//0.624;
1360   }
1361
1362   coefficients[0] = c1;
1363   coefficients[1] = c2;
1364
1365   Double_t dt = 0.1;
1366
1367   rates[0] = TMath::Exp(-dt/(r1));
1368   rates[1] = TMath::Exp(-dt/(r2));
1369
1370   Int_t i = 0;
1371   Int_t k = 0;
1372
1373   Double_t reminder[2];
1374   Double_t correction = 0.0;
1375   Double_t result     = 0.0;
1376
1377   // Attention: computation order is important
1378   for (k = 0; k < nexp; k++) {
1379     reminder[k] = 0.0;
1380   }
1381
1382   for (i = 0; i < n; i++) {
1383
1384     result    = (source[i] - correction);    // No rescaling
1385     target[i] = result;
1386
1387     for (k = 0; k < nexp; k++) {
1388       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1389     }
1390
1391     correction = 0.0;
1392     for (k = 0; k < nexp; k++) {
1393       correction += reminder[k];
1394     }
1395
1396   }
1397
1398 }
1399
1400 //_____________________________________________________________________________
1401 void AliTRDclusterizer::ResetRecPoints()
1402 {
1403   //
1404   // Resets the list of rec points
1405   //
1406
1407   if (fRecPoints) {
1408     fRecPoints->Delete();
1409     delete fRecPoints;
1410   }
1411 }
1412
1413 //_____________________________________________________________________________
1414 TClonesArray *AliTRDclusterizer::RecPoints()
1415 {
1416   //
1417   // Returns the list of rec points
1418   //
1419
1420   if (!fRecPoints) {
1421     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1422       // determine number of clusters which has to be allocated
1423       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1424       if(fReconstructor->IsHLT()) nclusters /= AliTRDgeometry::kNsector;
1425
1426       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1427     }
1428     //SetClustersOwner(kTRUE);
1429     AliTRDReconstructor::SetClusters(0x0);
1430   }
1431   return fRecPoints;
1432
1433 }
1434
1435 //_____________________________________________________________________________
1436 void AliTRDclusterizer::FillLUT()
1437 {
1438   //
1439   // Create the LUT
1440   //
1441
1442   const Int_t kNlut = 128;
1443
1444   fLUTbin = AliTRDgeometry::kNlayer * kNlut;
1445
1446   // The lookup table from Bogdan
1447   Float_t lut[AliTRDgeometry::kNlayer][kNlut] = {
1448     {
1449       0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1450       0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1451       0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1452       0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1453       0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1454       0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1455       0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1456       0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1457       0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1458       0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1459       0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1460       0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1461       0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1462       0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1463       0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1464       0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1465     },
1466     {
1467       0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1468       0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1469       0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1470       0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1471       0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1472       0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1473       0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1474       0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1475       0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1476       0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1477       0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1478       0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1479       0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1480       0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1481       0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1482       0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1483     },
1484     {
1485       0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1486       0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1487       0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1488       0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1489       0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1490       0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1491       0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1492       0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1493       0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1494       0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1495       0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1496       0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1497       0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1498       0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1499       0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1500       0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1501     },
1502     {
1503       0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1504       0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1505       0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1506       0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1507       0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1508       0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1509       0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1510       0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1511       0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1512       0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1513       0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1514       0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1515       0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1516       0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1517       0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1518       0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1519     },
1520     {
1521       0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1522       0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1523       0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1524       0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1525       0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1526       0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1527       0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1528       0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1529       0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1530       0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1531       0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1532       0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1533       0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1534       0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1535       0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1536       0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1537     },
1538     {
1539       0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1540       0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1541       0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1542       0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1543       0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1544       0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1545       0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1546       0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1547       0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1548       0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1549       0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1550       0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1551       0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1552       0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1553       0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1554       0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1555     }
1556   };
1557
1558   if (fLUT) {
1559     delete [] fLUT;
1560   }
1561   fLUT = new Double_t[fLUTbin];
1562
1563   for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1564     for (Int_t ilut  = 0; ilut  <  kNlut; ilut++  ) {
1565       fLUT[ilayer*kNlut+ilut] = lut[ilayer][ilut];
1566     }
1567   }
1568
1569 }
1570
1571 //_____________________________________________________________________________
1572 Double_t AliTRDclusterizer::LUTposition(Int_t ilayer
1573                                       , Double_t ampL
1574                                       , Double_t ampC
1575                                       , Double_t ampR) const
1576 {
1577   //
1578   // Calculates the cluster position using the lookup table.
1579   // Method provided by Bogdan Vulpescu.
1580   //
1581
1582   const Int_t kNlut = 128;
1583
1584   Double_t pos;
1585   Double_t x    = 0.0;
1586   Double_t xmin;
1587   Double_t xmax;
1588   Double_t xwid;
1589
1590   Int_t    side = 0;
1591   Int_t    ix;
1592
1593   Double_t xMin[AliTRDgeometry::kNlayer] = { 0.006492, 0.006377, 0.006258
1594                                           , 0.006144, 0.006030, 0.005980 };
1595   Double_t xMax[AliTRDgeometry::kNlayer] = { 0.960351, 0.965870, 0.970445
1596                                           , 0.974352, 0.977667, 0.996101 };
1597
1598   if      (ampL > ampR) {
1599     x    = (ampL - ampR) / ampC;
1600     side = -1;
1601   }
1602   else if (ampL < ampR) {
1603     x    = (ampR - ampL) / ampC;
1604     side = +1;
1605   }
1606
1607   if (ampL != ampR) {
1608
1609     xmin = xMin[ilayer] + 0.000005;
1610     xmax = xMax[ilayer] - 0.000005;
1611     xwid = (xmax - xmin) / 127.0;
1612
1613     if      (x < xmin) {
1614       pos = 0.0000;
1615     }
1616     else if (x > xmax) {
1617       pos = side * 0.5000;
1618     }
1619     else {
1620       ix  = (Int_t) ((x - xmin) / xwid);
1621       pos = side * fLUT[ilayer*kNlut+ix];
1622     }
1623
1624   }
1625   else {
1626
1627     pos = 0.0;
1628
1629   }
1630
1631   return pos;
1632
1633 }