TRD/AliTRDclusterizerV2.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id$ */
  18
  19 ///////////////////////////////////////////////////////////////////////////////
  20 //                                                                           //
  21 // TRD cluster finder                                                        //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include <TF1.h>
  26 #include <TTree.h>
  27 #include <TH1.h>
  28 #include <TFile.h>
  29
  30 #include "AliRunLoader.h"
  31 #include "AliLoader.h"
  32 #include "AliRawReader.h"
  33 #include "AliLog.h"
  34 #include "AliAlignObj.h"
  35
  36 #include "AliTRDclusterizerV2.h"
  37 #include "AliTRDgeometry.h"
  38 #include "AliTRDdataArrayF.h"
  39 #include "AliTRDdataArrayI.h"
  40 #include "AliTRDdigitsManager.h"
  41 #include "AliTRDrawData.h"
  42 #include "AliTRDcalibDB.h"
  43 #include "AliTRDSimParam.h"
  44 #include "AliTRDRecParam.h"
  45 #include "AliTRDcluster.h"
  46 #include "AliTRDtransform.h"
  47 #include "AliTRDSignalIndex.h"
  48 #include "AliTRDRawStream.h"
  49 #include "AliTRDRawStreamV2.h"
  50 #include "AliTRDfeeParam.h"
  51
  52 #include "Cal/AliTRDCalROC.h"
  53 #include "Cal/AliTRDCalDet.h"
  54
  55 ClassImp(AliTRDclusterizerV2)
  56
  57 //_____________________________________________________________________________
  58 AliTRDclusterizerV2::AliTRDclusterizerV2()
  59   :AliTRDclusterizer()
  60   ,fDigitsManager(NULL)
  61   ,fAddLabels(kTRUE)
  62   ,fRawVersion(2)
  63   ,fIndexesOut(NULL)
  64   ,fIndexesMaxima(NULL)
  65   ,fTransform(NULL)
  66 {
  67   //
  68   // AliTRDclusterizerV2 default constructor
  69   //
  70
  71   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
  72
  73 }
  74
  75 //_____________________________________________________________________________
  76 AliTRDclusterizerV2::AliTRDclusterizerV2(const Text_t *name, const Text_t *title)
  77   :AliTRDclusterizer(name,title)
  78   ,fDigitsManager(new AliTRDdigitsManager())
  79   ,fAddLabels(kTRUE)
  80   ,fRawVersion(2)
  81   ,fIndexesOut(NULL)
  82   ,fIndexesMaxima(NULL)
  83   ,fTransform(new AliTRDtransform(0))
  84 {
  85   //
  86   // AliTRDclusterizerV2 constructor
  87   //
  88
  89   fDigitsManager->CreateArrays();
  90
  91   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
  92
  93 }
  94
  95 //_____________________________________________________________________________
  96 AliTRDclusterizerV2::AliTRDclusterizerV2(const AliTRDclusterizerV2 &c)
  97   :AliTRDclusterizer(c)
  98   ,fDigitsManager(NULL)
  99   ,fAddLabels(kTRUE)
 100   ,fRawVersion(2)
 101   ,fIndexesOut(NULL)
 102   ,fIndexesMaxima(NULL)
 103   ,fTransform(NULL)
 104 {
 105   //
 106   // AliTRDclusterizerV2 copy constructor
 107   //
 108
 109 }
 110
 111 //_____________________________________________________________________________
 112 AliTRDclusterizerV2::~AliTRDclusterizerV2()
 113 {
 114   //
 115   // AliTRDclusterizerV2 destructor
 116   //
 117
 118   if (fDigitsManager)
 119     {
 120       delete fDigitsManager;
 121       fDigitsManager = NULL;
 122     }
 123
 124   if (fIndexesOut)
 125     {
 126       delete fIndexesOut;
 127       fIndexesOut    = NULL;
 128     }
 129
 130   if (fIndexesMaxima)
 131     {
 132       delete fIndexesMaxima;
 133       fIndexesMaxima = NULL;
 134     }
 135
 136   if (fTransform)
 137     {
 138       delete fTransform;
 139       fTransform     = NULL;
 140     }
 141
 142 }
 143
 144 //_____________________________________________________________________________
 145 AliTRDclusterizerV2 &AliTRDclusterizerV2::operator=(const AliTRDclusterizerV2 &c)
 146 {
 147   //
 148   // Assignment operator
 149   //
 150
 151   if (this != &c) ((AliTRDclusterizerV2 &) c).Copy(*this);
 152   return *this;
 153
 154 }
 155
 156 //_____________________________________________________________________________
 157 void AliTRDclusterizerV2::Copy(TObject &c) const
 158 {
 159   //
 160   // Copy function
 161   //
 162
 163   ((AliTRDclusterizerV2 &) c).fDigitsManager = NULL;
 164   ((AliTRDclusterizerV2 &) c).fAddLabels     = fAddLabels;
 165   ((AliTRDclusterizerV2 &) c).fRawVersion    = fRawVersion;
 166   ((AliTRDclusterizerV2 &) c).fIndexesOut    = NULL;
 167   ((AliTRDclusterizerV2 &) c).fIndexesMaxima = NULL;
 168   ((AliTRDclusterizerV2 &) c).fTransform     = NULL;
 169
 170   AliTRDclusterizer::Copy(c);
 171
 172 }
 173
 174 //_____________________________________________________________________________
 175 void AliTRDclusterizerV2::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
 176 {
 177   //
 178   // Reset the helper indexes
 179   //
 180
 181   if (fIndexesOut)
 182     {
 183       // carefull here - we assume that only row number may change - most probable
 184       if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
 185         fIndexesOut->ResetContent();
 186       else
 187         fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
 188                                            , indexesIn->GetNcol()
 189                                            , indexesIn->GetNtime());
 190     }
 191   else
 192     {
 193       fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
 194                                         , indexesIn->GetNcol()
 195                                         , indexesIn->GetNtime());
 196     }
 197
 198   if (fIndexesMaxima)
 199     {
 200       // carefull here - we assume that only row number may change - most probable
 201       if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
 202         {
 203           fIndexesMaxima->ResetContent();
 204         }
 205       else
 206         {
 207           fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
 208                                                 , indexesIn->GetNcol()
 209                                                 , indexesIn->GetNtime());
 210         }
 211     }
 212   else
 213     {
 214       fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
 215                                            , indexesIn->GetNcol()
 216                                            , indexesIn->GetNtime());
 217     }
 218
 219 }
 220
 221 //_____________________________________________________________________________
 222 Bool_t AliTRDclusterizerV2::ReadDigits()
 223 {
 224   //
 225   // Reads the digits arrays from the input aliroot file
 226   //
 227
 228   if (!fRunLoader) {
 229     AliError("No run loader available");
 230     return kFALSE;
 231   }
 232
 233   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
 234   if (!loader->TreeD()) {
 235     loader->LoadDigits();
 236   }
 237
 238   // Read in the digit arrays
 239   return (fDigitsManager->ReadDigits(loader->TreeD()));
 240
 241 }
 242
 243 //_____________________________________________________________________________
 244 Bool_t AliTRDclusterizerV2::ReadDigits(TTree *digitsTree)
 245 {
 246   //
 247   // Reads the digits arrays from the input tree
 248   //
 249
 250   // Read in the digit arrays
 251   return (fDigitsManager->ReadDigits(digitsTree));
 252
 253 }
 254
 255 //_____________________________________________________________________________
 256 Bool_t AliTRDclusterizerV2::ReadDigits(AliRawReader *rawReader)
 257 {
 258   //
 259   // Reads the digits arrays from the ddl file
 260   //
 261
 262   AliTRDrawData raw;
 263   fDigitsManager = raw.Raw2Digits(rawReader);
 264
 265   return kTRUE;
 266
 267 }
 268
 269 //_____________________________________________________________________________
 270 Bool_t AliTRDclusterizerV2::MakeClusters()
 271 {
 272   //
 273   // Creates clusters from digits
 274   //
 275
 276   // Propagate info from the digits manager
 277   if (fAddLabels == kTRUE)
 278     {
 279       fAddLabels = fDigitsManager->UsesDictionaries();
 280     }
 281
 282   Bool_t fReturn = kTRUE;
 283   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
 284     {
 285
 286       AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(i);
 287       // This is to take care of switched off super modules
 288       if (digitsIn->GetNtime() == 0)
 289         {
 290           continue;
 291         }
 292       digitsIn->Expand();
 293       AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
 294       if (indexes->IsAllocated() == kFALSE)
 295         {
 296           fDigitsManager->BuildIndexes(i);
 297         }
 298
 299       Bool_t fR = kFALSE;
 300       if (indexes->HasEntry())
 301         {
 302           if (fAddLabels)
 303             {
 304               for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
 305                 {
 306                   AliTRDdataArrayI *tracksIn = 0;
 307                   tracksIn = fDigitsManager->GetDictionary(i,iDict);
 308                   tracksIn->Expand();
 309                 }
 310             }
 311           fR = MakeClusters(i);
 312           fReturn = fR && fReturn;
 313         }
 314
 315       if (fR == kFALSE)
 316         {
 317           WriteClusters(i);
 318           ResetRecPoints();
 319         }
 320
 321       // No compress just remove
 322       fDigitsManager->RemoveDigits(i);
 323       fDigitsManager->RemoveDictionaries(i);
 324       fDigitsManager->ClearIndexes(i);
 325
 326     }
 327
 328   return fReturn;
 329
 330 }
 331
 332 //_____________________________________________________________________________
 333 Bool_t AliTRDclusterizerV2::Raw2Clusters(AliRawReader *rawReader)
 334 {
 335   //
 336   // Creates clusters from raw data
 337   //
 338
 339   AliTRDdataArrayI *digits = 0;
 340   AliTRDdataArrayI *track0 = 0;
 341   AliTRDdataArrayI *track1 = 0;
 342   AliTRDdataArrayI *track2 = 0;
 343
 344   AliTRDSignalIndex *indexes = 0;
 345
 346   // Create the digits manager
 347   if (!fDigitsManager)
 348     {
 349       fDigitsManager = new AliTRDdigitsManager();
 350       fDigitsManager->CreateArrays();
 351     }
 352
 353   AliTRDRawStreamV2 input(rawReader);
 354   input.SetRawVersion( fRawVersion );
 355   input.Init();
 356
 357   AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
 358
 359   // Loop through the digits
 360   Int_t lastdet = -1;
 361   Int_t det     =  0;
 362   Int_t it      =  0;
 363   while (input.Next())
 364     {
 365
 366       det = input.GetDet();
 367
 368       if (det != lastdet)
 369         {
 370
 371           if (lastdet != -1)
 372             {
 373               digits = fDigitsManager->GetDigits(lastdet);
 374               Bool_t iclusterBranch = kFALSE;
 375               if (indexes->HasEntry())
 376                 iclusterBranch = MakeClusters(lastdet);
 377               if (iclusterBranch == kFALSE)
 378                 {
 379                   WriteClusters(lastdet);
 380                   ResetRecPoints();
 381                 }
 382             }
 383
 384           if (digits)
 385             {
 386               fDigitsManager->RemoveDigits(lastdet);
 387               fDigitsManager->RemoveDictionaries(lastdet);
 388               fDigitsManager->ClearIndexes(lastdet);
 389             }
 390
 391           lastdet = det;
 392
 393           // Add a container for the digits of this detector
 394           digits = fDigitsManager->GetDigits(det);
 395           track0 = fDigitsManager->GetDictionary(det,0);
 396           track1 = fDigitsManager->GetDictionary(det,1);
 397           track2 = fDigitsManager->GetDictionary(det,2);
 398
 399           // Allocate memory space for the digits buffer
 400           if (digits->GetNtime() == 0)
 401             {
 402               //AliDebug(5, Form("Alloc digits for det %d", det));
 403               digits->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 404               track0->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 405               track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 406               track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
 407             }
 408
 409           indexes = fDigitsManager->GetIndexes(det);
 410           indexes->SetSM(input.GetSM());
 411           indexes->SetStack(input.GetStack());
 412           indexes->SetLayer(input.GetLayer());
 413           indexes->SetDetNumber(det);
 414           if (indexes->IsAllocated() == kFALSE)
 415             {
 416               indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins());
 417             }
 418
 419         }
 420
 421       for (it = 0; it < 3; it++)
 422         {
 423           if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() )
 424             {
 425               if (input.GetSignals()[it] > 0)
 426                 {
 427                   digits->SetDataUnchecked(input.GetRow(), input.GetCol(),
 428                                            input.GetTimeBin() + it, input.GetSignals()[it]);
 429
 430                   indexes->AddIndexTBin(input.GetRow(), input.GetCol(),
 431                                         input.GetTimeBin() + it);
 432                   track0->SetDataUnchecked(input.GetRow(), input.GetCol(),
 433                                            input.GetTimeBin() + it, 0);
 434                   track1->SetDataUnchecked(input.GetRow(), input.GetCol(),
 435                                            input.GetTimeBin() + it, 0);
 436                   track2->SetDataUnchecked(input.GetRow(), input.GetCol(),
 437                                            input.GetTimeBin() + it, 0);
 438                 }
 439             }
 440         }
 441
 442     }
 443
 444   if (lastdet != -1)
 445     {
 446       Bool_t iclusterBranch = kFALSE;
 447       if (indexes->HasEntry())
 448         {
 449           iclusterBranch = MakeClusters(lastdet);
 450         }
 451       if (iclusterBranch == kFALSE)
 452         {
 453           WriteClusters(lastdet);
 454           ResetRecPoints();
 455         }
 456       //MakeClusters(lastdet);
 457       if (digits)
 458         {
 459           fDigitsManager->RemoveDigits(lastdet);
 460           fDigitsManager->RemoveDictionaries(lastdet);
 461           fDigitsManager->ClearIndexes(lastdet);
 462         }
 463     }
 464
 465   delete fDigitsManager;
 466   fDigitsManager = NULL;
 467   return kTRUE;
 468
 469 }
 470
 471 //_____________________________________________________________________________
 472 Bool_t AliTRDclusterizerV2::Raw2ClustersChamber(AliRawReader *rawReader)
 473 {
 474   //
 475   // Creates clusters from raw data
 476   //
 477
 478   // Create the digits manager
 479   if (!fDigitsManager)
 480     {
 481       fDigitsManager = new AliTRDdigitsManager();
 482       fDigitsManager->CreateArrays();
 483     }
 484
 485   fDigitsManager->SetUseDictionaries(fAddLabels);
 486
 487   AliTRDRawStreamV2 input(rawReader);
 488   input.SetRawVersion( fRawVersion );
 489   input.Init();
 490
 491   AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
 492
 493   Int_t det    = 0;
 494   while ((det = input.NextChamber(fDigitsManager)) >= 0)
 495     {
 496       Bool_t iclusterBranch = kFALSE;
 497       if (fDigitsManager->GetIndexes(det)->HasEntry())
 498         {
 499           iclusterBranch = MakeClusters(det);
 500         }
 501       if (iclusterBranch == kFALSE)
 502         {
 503           WriteClusters(det);
 504           ResetRecPoints();
 505         }
 506       fDigitsManager->RemoveDigits(det);
 507       fDigitsManager->RemoveDictionaries(det);
 508       fDigitsManager->ClearIndexes(det);
 509     }
 510
 511   delete fDigitsManager;
 512   fDigitsManager = NULL;
 513   return kTRUE;
 514
 515 }
 516
 517 //_____________________________________________________________________________
 518 Bool_t AliTRDclusterizerV2::MakeClusters(Int_t det)
 519 {
 520   //
 521   // Generates the cluster.
 522   //
 523
 524   // Get the digits
 525   //   digits should be expanded beforehand!
 526   //   digitsIn->Expand();
 527   AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(det);
 528
 529   // This is to take care of switched off super modules
 530   if (digitsIn->GetNtime() == 0)
 531     {
 532       return kFALSE;
 533     }
 534
 535   AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
 536   if (indexesIn->IsAllocated() == kFALSE)
 537     {
 538       AliError("Indexes do not exist!");
 539       return kFALSE;
 540     }
 541
 542   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
 543   if (!calibration)
 544     {
 545       AliFatal("No AliTRDcalibDB instance available\n");
 546       return kFALSE;
 547     }
 548
 549   AliTRDRecParam *recParam    = AliTRDRecParam::Instance();
 550   if (!recParam)
 551     {
 552       AliError("No AliTRDRecParam instance available\n");
 553       return kFALSE;
 554     }
 555
 556   // ADC thresholds
 557   // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
 558   Float_t ADCthreshold   = 0;
 559
 560   // Threshold value for the maximum
 561   Float_t maxThresh      = recParam->GetClusMaxThresh();
 562   // Threshold value for the digit signal
 563   Float_t sigThresh      = recParam->GetClusSigThresh();
 564
 565   // Iteration limit for unfolding procedure
 566   const Float_t kEpsilon = 0.01;
 567   const Int_t   kNclus   = 3;
 568   const Int_t   kNsig    = 5;
 569
 570   Int_t    iUnfold       = 0;
 571   Double_t ratioLeft     = 1.0;
 572   Double_t ratioRight    = 1.0;
 573
 574   Double_t padSignal[kNsig];
 575   Double_t clusterSignal[kNclus];
 576
 577   Int_t icham = indexesIn->GetChamber();
 578   Int_t iplan = indexesIn->GetPlane();
 579   Int_t isect = indexesIn->GetSM();
 580
 581   // Start clustering in the chamber
 582
 583   Int_t idet  = AliTRDgeometry::GetDetector(iplan,icham,isect);
 584   if (idet != det)
 585     {
 586       AliError("Strange Detector number Missmatch!");
 587       return kFALSE;
 588     }
 589
 590   // TRD space point transformation
 591   fTransform->SetDetector(det);
 592
 593   Int_t    ilayer  = AliGeomManager::kTRD1 + iplan;
 594   Int_t    imodule = icham + AliTRDgeometry::Ncham() * isect;
 595   UShort_t volid   = AliGeomManager::LayerToVolUID(ilayer,imodule);
 596
 597   Int_t nColMax    = digitsIn->GetNcol();
 598   Int_t nTimeTotal = digitsIn->GetNtime();
 599
 600   // Detector wise calibration object for the gain factors
 601   const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();
 602   // Calibration object with pad wise values for the gain factors
 603   AliTRDCalROC       *calGainFactorROC      = calibration->GetGainFactorROC(idet);
 604   // Calibration value for chamber wise gain factor
 605   Float_t             calGainFactorDetValue = calGainFactorDet->GetValue(idet);
 606
 607   Int_t nClusters = 0;
 608
 609   AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow()
 610                                                     ,digitsIn->GetNcol()
 611                                                     ,digitsIn->GetNtime());
 612
 613   ResetHelperIndexes(indexesIn);
 614
 615   // Apply the gain and the tail cancelation via digital filter
 616   TailCancelation(digitsIn
 617                  ,digitsOut
 618                  ,indexesIn
 619                  ,fIndexesOut
 620                  ,nTimeTotal
 621                  ,ADCthreshold
 622                  ,calGainFactorROC
 623                  ,calGainFactorDetValue);
 624
 625   Int_t row  = 0;
 626   Int_t col  = 0;
 627   Int_t time = 0;
 628   Int_t iPad = 0;
 629
 630   fIndexesOut->ResetCounters();
 631   while (fIndexesOut->NextRCTbinIndex(row, col, time))
 632     {
 633
 634       Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
 635
 636       // Look for the maximum
 637       if (signalM >= maxThresh)
 638         {
 639
 640           if (col + 1 >= nColMax || col-1 < 0)
 641             continue;
 642
 643           Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
 644           Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
 645
 646           if ((TMath::Abs(signalL) <= signalM) &&
 647               (TMath::Abs(signalR) <  signalM))
 648             {
 649               if ((TMath::Abs(signalL) >= sigThresh) ||
 650                   (TMath::Abs(signalR) >= sigThresh))
 651                 {
 652                   // Maximum found, mark the position by a negative signal
 653                   digitsOut->SetDataUnchecked(row,col,time,-signalM);
 654                   fIndexesMaxima->AddIndexTBin(row,col,time);
 655                 }
 656             }
 657
 658         }
 659
 660     }
 661
 662   // The index to the first cluster of a given ROC
 663   Int_t firstClusterROC = -1;
 664   // The number of cluster in a given ROC
 665   Int_t nClusterROC     =  0;
 666
 667   // Now check the maxima and calculate the cluster position
 668   fIndexesMaxima->ResetCounters();
 669   while (fIndexesMaxima->NextRCTbinIndex(row, col, time))
 670     {
 671
 672       // Maximum found ?
 673       if (digitsOut->GetDataUnchecked(row,col,time) < 0.0)
 674         {
 675
 676           for (iPad = 0; iPad < kNclus; iPad++)
 677             {
 678               Int_t iPadCol = col - 1 + iPad;
 679               clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
 680             }
 681
 682           // Count the number of pads in the cluster
 683           Int_t nPadCount = 0;
 684           Int_t ii;
 685           // Look to the left
 686           ii = 0;
 687           while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii  ,time)) >= sigThresh)
 688             {
 689               nPadCount++;
 690               ii++;
 691               if (col-ii   <        0) break;
 692             }
 693           // Look to the right
 694           ii = 0;
 695           while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh)
 696             {
 697               nPadCount++;
 698               ii++;
 699               if (col+ii+1 >= nColMax) break;
 700             }
 701           nClusters++;
 702
 703           // Look for 5 pad cluster with minimum in the middle
 704           Bool_t fivePadCluster = kFALSE;
 705           if (col < (nColMax - 3))
 706             {
 707               if (digitsOut->GetDataUnchecked(row,col+2,time) < 0)
 708                 {
 709                   fivePadCluster = kTRUE;
 710                 }
 711               if ((fivePadCluster) && (col < (nColMax - 5)))
 712                 {
 713                   if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh)
 714                     {
 715                       fivePadCluster = kFALSE;
 716                     }
 717                 }
 718               if ((fivePadCluster) && (col >             1))
 719                 {
 720                   if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh)
 721                     {
 722                       fivePadCluster = kFALSE;
 723                     }
 724                 }
 725             }
 726
 727           // 5 pad cluster
 728           // Modify the signal of the overlapping pad for the left part
 729           // of the cluster which remains from a previous unfolding
 730           if (iUnfold)
 731             {
 732               clusterSignal[0] *= ratioLeft;
 733               iUnfold = 0;
 734             }
 735
 736           // Unfold the 5 pad cluster
 737           if (fivePadCluster)
 738             {
 739               for (iPad = 0; iPad < kNsig; iPad++)
 740                 {
 741                   padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
 742                                                                           ,col-1+iPad
 743                                                                           ,time));
 744                 }
 745               // Unfold the two maxima and set the signal on
 746               // the overlapping pad to the ratio
 747               ratioRight        = Unfold(kEpsilon,iplan,padSignal);
 748               ratioLeft         = 1.0 - ratioRight;
 749               clusterSignal[2] *= ratioRight;
 750               iUnfold = 1;
 751             }
 752
 753           // The position of the cluster in COL direction relative to the center pad (pad units)
 754           Double_t clusterPosCol = 0.0;
 755           if (recParam->LUTOn())
 756             {
 757               // Calculate the position of the cluster by using the
 758               // lookup table method
 759               clusterPosCol = recParam->LUTposition(iplan
 760                                                    ,clusterSignal[0]
 761                                                    ,clusterSignal[1]
 762                                                    ,clusterSignal[2]);
 763             }
 764           else
 765             {
 766               // Calculate the position of the cluster by using the
 767               // center of gravity method
 768               for (Int_t i = 0; i < kNsig; i++)
 769                 {
 770                   padSignal[i] = 0.0;
 771                 }
 772               padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col  ,time)); // Central pad
 773               padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left    pad
 774               padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right   pad
 775               if ((col >           2) &&
 776                   (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1]))
 777                 {
 778                   padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
 779                 }
 780               if ((col < nColMax - 3) &&
 781                   (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3]))
 782                 {
 783                   padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
 784                 }
 785               clusterPosCol = GetCOG(padSignal);
 786             }
 787
 788           // Store the amplitudes of the pads in the cluster for later analysis
 789           Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
 790           for (Int_t jPad = col-3; jPad <= col+3; jPad++)
 791             {
 792               if ((jPad <          0) ||
 793                   (jPad >= nColMax-1))
 794                 {
 795                   continue;
 796                 }
 797               signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
 798             }
 799
 800           // Transform the local cluster coordinates into calibrated
 801           // space point positions defined in the local tracking system.
 802           // Here the calibration for T0, Vdrift and ExB is applied as well.
 803           Double_t clusterXYZ[6];
 804           clusterXYZ[0] = clusterPosCol;
 805           clusterXYZ[1] = clusterSignal[0];
 806           clusterXYZ[2] = clusterSignal[1];
 807           clusterXYZ[3] = clusterSignal[2];
 808           clusterXYZ[4] = 0.0;
 809           clusterXYZ[5] = 0.0;
 810           Int_t    clusterRCT[3];
 811           clusterRCT[0] = row;
 812           clusterRCT[1] = col;
 813           clusterRCT[2] = 0;
 814           fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),0);
 815
 816           // Add the cluster to the output array
 817           // The track indices will be stored later
 818           Float_t clusterPos[3];
 819           clusterPos[0] = clusterXYZ[0];
 820           clusterPos[1] = clusterXYZ[1];
 821           clusterPos[2] = clusterXYZ[2];
 822           Float_t clusterSig[2];
 823           clusterSig[0] = clusterXYZ[4];
 824           clusterSig[1] = clusterXYZ[5];
 825           Double_t clusterCharge  = clusterXYZ[3];
 826           Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
 827           AliTRDcluster *cluster = new AliTRDcluster(idet
 828                                                     ,clusterCharge
 829                                                     ,clusterPos
 830                                                     ,clusterSig
 831                                                     ,0x0
 832                                                     ,((Char_t) nPadCount)
 833                                                     ,signals
 834                                                     ,((UChar_t) col)
 835                                                     ,((UChar_t) row)
 836                                                     ,((UChar_t) time)
 837                                                     ,clusterTimeBin
 838                                                     ,clusterPosCol
 839                                                     ,volid);
 840
 841           // Temporarily store the row, column and time bin of the center pad
 842           // Used to later on assign the track indices
 843           cluster->SetLabel( row,0);
 844           cluster->SetLabel( col,1);
 845           cluster->SetLabel(time,2);
 846
 847           RecPoints()->Add(cluster);
 848
 849           // Store the index of the first cluster in the current ROC
 850           if (firstClusterROC < 0)
 851             {
 852               firstClusterROC = RecPoints()->GetEntriesFast() - 1;
 853             }
 854
 855           // Count the number of cluster in the current ROC
 856           nClusterROC++;
 857
 858         } // if: Maximum found ?
 859
 860     }
 861
 862   delete digitsOut;
 863
 864   if (fAddLabels)
 865     {
 866       AddLabels(idet, firstClusterROC, nClusterROC);
 867     }
 868
 869   // Write the cluster and reset the array
 870   WriteClusters(idet);
 871   ResetRecPoints();
 872
 873   return kTRUE;
 874
 875 }
 876
 877 //_____________________________________________________________________________
 878 Bool_t AliTRDclusterizerV2::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
 879 {
 880   //
 881   // Add the track indices to the found clusters
 882   //
 883
 884   const Int_t   kNclus  = 3;
 885   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
 886   const Int_t   kNtrack = kNdict * kNclus;
 887
 888   Int_t iClusterROC = 0;
 889
 890   Int_t row  = 0;
 891   Int_t col  = 0;
 892   Int_t time = 0;
 893   Int_t iPad = 0;
 894
 895   // Temporary array to collect the track indices
 896   Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
 897
 898   // Loop through the dictionary arrays one-by-one
 899   // to keep memory consumption low
 900   AliTRDdataArrayI *tracksIn = 0;
 901   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
 902
 903     // tracksIn should be expanded beforehand!
 904     tracksIn = fDigitsManager->GetDictionary(idet,iDict);
 905
 906     // Loop though the clusters found in this ROC
 907     for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
 908
 909       AliTRDcluster *cluster = (AliTRDcluster *)
 910         RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
 911       row  = cluster->GetLabel(0);
 912       col  = cluster->GetLabel(1);
 913       time = cluster->GetLabel(2);
 914
 915       for (iPad = 0; iPad < kNclus; iPad++) {
 916         Int_t iPadCol = col - 1 + iPad;
 917         Int_t index   = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
 918         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
 919       }
 920
 921     }
 922
 923   }
 924
 925   // Copy the track indices into the cluster
 926   // Loop though the clusters found in this ROC
 927   for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
 928
 929     AliTRDcluster *cluster = (AliTRDcluster *)
 930       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
 931     cluster->SetLabel(-9999,0);
 932     cluster->SetLabel(-9999,1);
 933     cluster->SetLabel(-9999,2);
 934
 935     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
 936
 937   }
 938
 939   delete [] idxTracks;
 940
 941   return kTRUE;
 942
 943 }
 944
 945 //_____________________________________________________________________________
 946 Double_t AliTRDclusterizerV2::GetCOG(Double_t signal[5])
 947 {
 948   //
 949   // Get COG position
 950   // Used for clusters with more than 3 pads - where LUT not applicable
 951   //
 952
 953   Double_t sum = signal[0]
 954                + signal[1]
 955                + signal[2]
 956                + signal[3]
 957                + signal[4];
 958
 959   Double_t res = (0.0 * (-signal[0] + signal[4])
 960                       + (-signal[1] + signal[3])) / sum;
 961
 962   return res;
 963
 964 }
 965
 966 //_____________________________________________________________________________
 967 Double_t AliTRDclusterizerV2::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
 968 {
 969   //
 970   // Method to unfold neighbouring maxima.
 971   // The charge ratio on the overlapping pad is calculated
 972   // until there is no more change within the range given by eps.
 973   // The resulting ratio is then returned to the calling method.
 974   //
 975
 976   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 977   if (!calibration) {
 978     AliError("No AliTRDcalibDB instance available\n");
 979     return kFALSE;
 980   }
 981
 982   Int_t   irc                = 0;
 983   Int_t   itStep             = 0;                 // Count iteration steps
 984
 985   Double_t ratio             = 0.5;               // Start value for ratio
 986   Double_t prevRatio         = 0.0;               // Store previous ratio
 987
 988   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
 989   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
 990   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
 991
 992   // Start the iteration
 993   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
 994
 995     itStep++;
 996     prevRatio = ratio;
 997
 998     // Cluster position according to charge ratio
 999     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
1000                       / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1001     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1002                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1003
1004     // Set cluster charge ratio
1005     irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
1006     Double_t ampLeft  = padSignal[1] / newSignal[1];
1007     irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
1008     Double_t ampRight = padSignal[3] / newSignal[1];
1009
1010     // Apply pad response to parameters
1011     irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
1012     irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
1013
1014     // Calculate new overlapping ratio
1015     ratio = TMath::Min((Double_t) 1.0
1016                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1017
1018   }
1019
1020   return ratio;
1021
1022 }
1023
1024 //_____________________________________________________________________________
1025 void AliTRDclusterizerV2::TailCancelation(AliTRDdataArrayI *digitsIn
1026                                         , AliTRDdataArrayF *digitsOut
1027                                         , AliTRDSignalIndex *indexesIn
1028                                         , AliTRDSignalIndex *indexesOut
1029                                         , Int_t nTimeTotal
1030                                         , Float_t ADCthreshold
1031                                         , AliTRDCalROC *calGainFactorROC
1032                                         , Float_t calGainFactorDetValue)
1033 {
1034   //
1035   // Applies the tail cancelation and gain factors:
1036   // Transform digitsIn to digitsOut
1037   //
1038
1039   Int_t iRow  = 0;
1040   Int_t iCol  = 0;
1041   Int_t iTime = 0;
1042
1043   AliTRDRecParam *recParam = AliTRDRecParam::Instance();
1044   if (!recParam)
1045     {
1046       AliError("No AliTRDRecParam instance available\n");
1047       return;
1048     }
1049
1050   Double_t *inADC  = new Double_t[nTimeTotal];  // ADC data before tail cancellation
1051   Double_t *outADC = new Double_t[nTimeTotal];  // ADC data after tail cancellation
1052   indexesIn->ResetCounters();
1053   while (indexesIn->NextRCIndex(iRow, iCol))
1054     {
1055       Float_t  calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1056       Double_t gain                  = calGainFactorDetValue
1057                                      * calGainFactorROCValue;
1058
1059       for (iTime = 0; iTime < nTimeTotal; iTime++)
1060         {
1061           // Apply gain gain factor
1062           inADC[iTime]   = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1063           inADC[iTime]  /= gain;
1064           outADC[iTime]  = inADC[iTime];
1065         }
1066
1067       // Apply the tail cancelation via the digital filter
1068       if (recParam->TCOn())
1069         {
1070           DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp());
1071         }
1072
1073       indexesIn->ResetTbinCounter();
1074       while (indexesIn->NextTbinIndex(iTime))
1075         {
1076           // Store the amplitude of the digit if above threshold
1077           if (outADC[iTime] > ADCthreshold)
1078             {
1079               digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1080               indexesOut->AddIndexTBin(iRow,iCol,iTime);
1081             }
1082         } //while itime
1083
1084     }//while irow icol
1085
1086   delete [] inADC;
1087   delete [] outADC;
1088
1089   return;
1090
1091 }
1092
1093 //_____________________________________________________________________________
1094 void AliTRDclusterizerV2::DeConvExp(Double_t *source, Double_t *target
1095                                   , Int_t n, Int_t nexp)
1096 {
1097   //
1098   // Tail cancellation by deconvolution for PASA v4 TRF
1099   //
1100
1101   Double_t rates[2];
1102   Double_t coefficients[2];
1103
1104   // Initialization (coefficient = alpha, rates = lambda)
1105   Double_t r1 = 1.0;
1106   Double_t r2 = 1.0;
1107   Double_t c1 = 0.5;
1108   Double_t c2 = 0.5;
1109
1110   if (nexp == 1) {   // 1 Exponentials
1111     r1 = 1.156;
1112     r2 = 0.130;
1113     c1 = 0.066;
1114     c2 = 0.000;
1115   }
1116   if (nexp == 2) {   // 2 Exponentials
1117     r1 = 1.156;
1118     r2 = 0.130;
1119     c1 = 0.114;
1120     c2 = 0.624;
1121   }
1122
1123   coefficients[0] = c1;
1124   coefficients[1] = c2;
1125
1126   Double_t dt = 0.1;
1127
1128   rates[0] = TMath::Exp(-dt/(r1));
1129   rates[1] = TMath::Exp(-dt/(r2));
1130
1131   Int_t i = 0;
1132   Int_t k = 0;
1133
1134   Double_t reminder[2];
1135   Double_t correction = 0.0;
1136   Double_t result     = 0.0;
1137
1138   // Attention: computation order is important
1139   for (k = 0; k < nexp; k++) {
1140     reminder[k] = 0.0;
1141   }
1142
1143   for (i = 0; i < n; i++) {
1144
1145     result    = (source[i] - correction);    // No rescaling
1146     target[i] = result;
1147
1148     for (k = 0; k < nexp; k++) {
1149       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1150     }
1151
1152     correction = 0.0;
1153     for (k = 0; k < nexp; k++) {
1154       correction += reminder[k];
1155     }
1156
1157   }
1158
1159 }