TRD/AliTRDclusterizer.cxx

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