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